[SciPy-user] += Operator and Slicing of Arrays

Emmanuelle Gouillart emmanuelle.gouillart@normalesup....
Fri Jun 12 13:28:13 CDT 2009


Hi Josh,

what kind of problem do you have exactly ? Do you have trouble
implementing the computation you describe, or do you get unexpected
results?

If I understood well what you want to do, you cannot use directly use +=
with fancy indexing (quantity[tri_idx_plus,...]) because the repeated
elements will be incremented just once (see
http://www.scipy.org/Tentative_NumPy_Tutorial#head-0dffc419afa7d77d51062d40d2d84143db8216c2
for more details).

However, I think you can solve your problem by using a weighted
histogram. Using your notations 

weights = edge_unknown[:len(tri_idx_plus),...] *  
	basis_p[:len(tri_idx_plus),...]
histogram_values = np.histogram(tri_idx_plus,
	np.arange(tri_idx_plus.max() +2), weights=weights)
unique_plus = np.unique1d(tri_idx_plus)
quantity[unique_plus,...] = histogram_values[0][unique_plus]

The weighted histogram allows you to make the sums corresponding to each
triangle.

Here is an example

>>> quantity = np.zeros(10)
>>> tri_idx_plus = np.array([0, 3, 6, 3, 2, 4, 1, 4])
>>> weights = 2*tri_idx_plus + 1
>>> histogram_values = np.histogram(tri_idx_plus,
        np.arange(tri_idx_plus.max() +2), weights=weights)
>>> unique_plus = np.unique1d(tri_idx_plus)
>>> quantity[unique_plus,...] = histogram_values[0][unique_plus]

Actually, this may only work with positive values of weights (not
checked)...

Please tell us if this meets your needs or not.

Cheers,

Emmanuelle

On Fri, Jun 12, 2009 at 01:14:16PM -0400, Josh Lawrence wrote:
> Greetings,

> I am in need of some help. I am trying to use the += operator to sum  
> over edge elements on a triangular mesh. Each edge has an unknown  
> associated with it. After solving for the unknowns, I am trying to  
> compute a quantity at the centroid of all triangles in the mesh. On a  
> flat surface, each edge will be connected to either one or two  
> triangles. The orientation of the edge and the normals of the  
> triangles determines whether each triangle attached to the edge is a  
> "plus" or "minus" triangle for that edge. It is possible for one  
> triangle to be referenced three times as a "plus" triangle, three  
> times as a "minus" triangle or any combination of "plus" and  
> "minus" (1 and 2 or 2 and 1, respectively).

> I have a variable tri_idx which relates the edges to the "plus" and  
> "minus" triangles. I then compute the quantity at the centroid for the  
> "plus" triangle and "minus" triangle attached to each edge. An example  
> follows:

> tri_idx_plus = [0 3 6 3 2 4 1 4]
> tri_idx_minus = [1 2 5 3 6 0 1 4]

> quantity[tri_idx_plus,...] += edge_unknown[:len(tri_idx_plus),...] *  
> basis_p[:len(tri_idx_plus),...]
> quantity[tri_idx_minus,...] += edge_unknown[:len(tri_idx_minus)...] *  
> basis_m[:len(tri_idx_minus),...]

> where basis_p and basis_m are basis functions that expand the unknown  
> of each edge into a surface function over the "plus" or "minus"  
> triangle.

> I am pretty sure the problem I am encountering is that tri_idx_plus  
> mentions indices 3 and 4 twice and tri_idx_minus contains index 1  
> twice. Is there a way of doing this operation without reverting to  
> looping over each edge (read: not doing this the slow way).

> Thanks in advance!

> Josh Lawrence
> Ph.D. Student
> Clemson University

> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user


More information about the SciPy-user mailing list