Triangle project helper (Cell to point data converter)

[ZAARAOUI999]

Hello, Mr. Andreas, I think that the problem of the missing triangle project helper can be solved by this, take a look.

import felupe as fem
import numpy as np

import pygmsh


h = 1.0
L = 1.0

with pygmsh.geo.Geometry() as geom:
    
    rectangle = geom.add_rectangle(0,h,0,L, 0, mesh_size=h/10)
    mesh_ = geom.generate_mesh(dim=2, algorithm=2)
    

m = fem.Mesh(mesh_.points[:,:2], mesh_.cells[1].data, mesh_.cells[1].type)

region = fem.RegionTriangle(m)

displacement = fem.Field(region, dim=2)

field = fem.FieldContainer([displacement])

umat = fem.NeoHooke(mu=1, bulk=1)

solid = fem.SolidBody(umat, field)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)

move = fem.math.linsteps([0,.2],num=10)

ramp = {boundaries['move']:move}

step = fem.Step([solid], ramp, boundaries)

job = fem.Job([step])
job.evaluate(verbose=False) 


def tri_project(values, region):
    
    mesh = region.mesh.copy()
    dim = int(np.product(values.shape[:-2]))
    u = values.T.reshape(-1, dim)
    u = u.mean(-1)
    cells = mesh.cells
    points_data = np.zeros(mesh.npoints)
    cells_ids = np.arange(cells.shape[0])
    
    for i in range(mesh.npoints):
        cd = np.isin(cells, cells[cells==i])
        point_cells = cells_ids[np.any(cd,axis=-1)]
        point_data = u.take(point_cells).mean(-1)
        points_data[i] = point_data
        
    return points_data

# displacements along the X axis    
U1 = solid.field.extract(grad=False)[0][0]
points_data = tri_project(U1, region)

Then, I got this

I hope that this can help.

[ZAARAOUI999]

This is worked also with tetra mesh

[adtzlr]

Thank you very much! I'll have a look at it during the weekend.

[adtzlr]

Btw, which plotting library are you using?

[ZAARAOUI999]

I'm using PyVista, it's a super visualization tool.

[ZAARAOUI999]

I'm working on a visualization GUI tool for FElupe. So, the user may be forced to work with tetra or tri meshes that's why I had a need to the tri\tetra project helper. Also, I hope that an interaction module will be added to FElupe in the next days 😁 , this will be amazing!

[adtzlr]

I enhanced the code slightly to support tensor-valued data (like the original project-helper does also). However, I'd like to get rid of the for-loop but I have no idea for that right now.

def tri_project(values, region):
    """Projection of scalar or vectorial values taken from the mean of the 
    quadrature point data of cells to mesh-points (for triangles and tetrahedrons).
    """
    
    mesh = region.mesh
    cells = mesh.cells
    
    # 1d-reshaped values at cells
    dim = int(np.product(values.shape[:-2]))
    
    # mean of values at quadrature points of cells
    # reshape mean-values of cells to (ncells, dim)
    u = values.T.mean(1).reshape(-1, dim)

    # init values at mesh-points
    values_points = np.zeros((mesh.npoints, dim))
    
    # loop over mesh-points
    for i, point_data in enumerate(values_points):
        
        # mask cells which contain the given point
        mask = np.any(np.isin(cells, i), axis=1)
        
        # collect all cells which are connected to the given point
        point_cells = np.arange(mesh.ncells)[mask]
        
        # combine connected cell values to mean value at the given point
        values_points[i] = u[point_cells].mean(0)
        
    return values_points

[adtzlr]

This function has the potential to work will all kind of cells (mean values per cell of all values located at quadrature points). For quadrature rules with non uniform weights, the weighted mean must be evaluated. Probably we should call it project_cells_mean() or mean_project()?

[ZAARAOUI999]

I'm a lazy programmer so I will choose the second one, but at the end, it's your choice Mr. Andreas, just pick a one and let FElupe move another step 👍

[ZAARAOUI999]

I fired the for-loop. But there is no enhancement at the level of execution time, take a look;

A part of the code:

def tri_project(values, region):
    """Projection of scalar or vectorial values taken from the mean of the 
    quadrature point data of cells to mesh-points (for triangles and tetrahedrons).
    """
    
    mesh = region.mesh
    cells = mesh.cells
    
    # 1d-reshaped values at cells
    dim = int(np.product(values.shape[:-2]))
    
    # mean of values at quadrature points of cells
    # reshape mean-values of cells to (ncells, dim)
    u = values.T.mean(1).reshape(-1, dim)

    # init values at mesh-points
    values_points = np.zeros((mesh.npoints, dim))
    
    # loop over mesh-points
    for i, point_data in enumerate(values_points):
        
        # mask cells which contain the given point
        mask = np.any(np.isin(cells, i), axis=1)
        
        # collect all cells which are connected to the given point
        point_cells = np.arange(mesh.ncells)[mask]
        
        # combine connected cell values to mean value at the given point
        values_points[i] = u[point_cells].mean(0)
        
    return values_points

def tri_project_modified(values, region):
    """..."""
    
    mesh = region.mesh
    cells = mesh.cells
    
    # 1d-reshaped values at cells
    dim = int(np.product(values.shape[:-2]))    
    # mean of values at quadrature points of cells
    # reshape mean-values of cells to (ncells, dim)
    u = values.T.mean(1).reshape(-1, dim)
    
    points = np.arange(mesh.npoints)
    
    def compute_mean_value(point):
        # mask cells which contain the given point
        mask = np.any(np.isin(cells, point), axis=1)
        
        # mask cells which contain the given point
        point_cells = np.arange(mesh.ncells)[mask]
        
        # combine connected cell values to mean value at the given point
        value = u[point_cells].mean(0)
        
        return value
    
    # get mean values
    values = list(map(compute_mean_value, points))
    
    return values


def f_time(f,*args):
    t0 = time.time()
    result = f(*args)
    tf = time.time()
    t = tf - t0
    return result, t

# displacements along the X axis    
U1 = solid.field.extract(grad=False)[0][0]
try: 
    print(str(region.mesh.npoints)+' points')
    points_data_, t1 = f_time(tri_project, U1, region)
    points_data, t2 = f_time(tri_project_modified, U1, region)
    print('unmodified function: %.6f s'%t1)
    print('modified function: %.6f s'%t2)
except Exception as e:
    print(e)

then I got;

3435 points
unmodified function: 1.421936 s
modified function: 1.419627 s

It stills no meaningful difference.

[adtzlr]

Thank's for your effort! It's because of map, list-comprehensions and for-loops are all more or less equal in execution times. I thought of some NumPy-magic.

FYI: FElupe relies a lot on vectorized operations, especially on np.einsum() with leading or trailing axes. This way you get a fantastic performance on applying the same operation on lots of elements again and again, even in interpreted languages like Python (without any compilation). Best of it, NumPy often takes care of parallel evaluations of underlying dot-products. A small example:

import numpy as np

F = np.random.rand(3, 3, 100000)
C = np.zeros_like(F)

# this is slow
for i in range(3):
    for j in range(3):
        for k in range(3):
            for a in range(100000):
                C[i, j, a] += F[k, i, a] * F[k, j, a]

transpose = lambda x: np.einsum("ij...->ji...", x)
dot = lambda x, y: np.einsum("ik...,kj...->ij...", x, y)

# much much faster
C2 = dot(transpose(F), F)

Of course this is not super-critical for the projection, but it would be nice to eliminate the for-loop.

[adtzlr]

@ZAARAOUI999 But your approach is not wasted time. Often this is very useful, e.g. when using something like joblib (though not used by FElupe) to evaluate a for-loop in parallel.

[ZAARAOUI999]

the joblib library looks like an interested lib, I will test it, thank you for all 🙂.

[adtzlr]

I created a project helper without a loop. First of all: it is not faster 🐢. Second, it is a memory beast 📚 - thanks to the huge array u * mask.

def tri_project_no_loop(values, region):
    """Projection of scalar or vectorial values taken from the mean of the 
    quadrature point data of cells to mesh-points (for triangles and tetrahedrons).
    """
    
    mesh = region.mesh
    
    # 1d-reshaped values at cells
    dim = int(np.product(values.shape[:-2]))
    
    # mean of values at quadrature points of cells
    u = values.mean(-2).reshape(dim, mesh.ncells, 1)
    
    # mask cells which contain a specific point of shape (ncells, npoints)
    mask = (mesh.cells[..., np.newaxis] == np.arange(mesh.npoints)).any(axis=1)

    # # obtain values on mesh-points
    return ((u * mask).sum(1) / mask.astype(bool).sum(0)).T

[adtzlr]

@ZAARAOUI999, please have quick look at the proposed implementation #329 if the solution is okay for you. Many thanks!

[ZAARAOUI999]

Sorry for being late, I can't find out where is the sparsematrix , I tried to set average to False but the result was an array that its dim didn't match with both of npoints and ncells. Can you give a quik exemple?

[ZAARAOUI999]

It's worked for average=True with a super speed (I found the sparsematrix).That's great 👍

[adtzlr]

Not averaging means the point results are not combined, the result operates on a disconnected mesh. No. of points is of course different because it is len(region.element.points) * mesh.ncells.

[ZAARAOUI999]

okey, that's clear now, thanks a lot for speeding up the helper