In this notebook, we illustrate the epithelium generation strategy.
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import graph_tool.all as gt
import pandas as pd
import seaborn as sns
import IPython.display as disp
import leg_joint
from leg_joint.epithelium.generation import base_grid, cylindrical
from leg_joint.epithelium.generation import (vertex_data, cell_data,
edge_data, junction_data,
face_data)
/home/guillaume/anaconda/envs/python3/lib/python3.4/site-packages/graph_tool/draw/cairo_draw.py:1318: RuntimeWarning: Error importing Gtk module: No module named 'gi'; GTK+ drawing will not work. warnings.warn(msg, RuntimeWarning)
epithelium.generation
documentation¶This module contains the elements necessary to build an epithelium.
vertex_data, cell_data, edge_data, junction_data and face_data are dictionnaries of the form:
{'data_name': (default_value, data_type)}.
They are used to instanciate all the necessary property maps over the graph.
vertex_data and edge_data contain data unspecific to the vertex or edge type (here cell vertices, junction edges, and triangular faces). Specific data for cells and junction edges will be instanciated as PropertyMaps of the corresponding GraphView, except for faces. Note that graphviews inherits the parent graphs propertymaps at instanciation time.
It is recommanded to add new properties by specifying them here, although they can also be specified when Epithelium is instanciated. Adding properties dynamically is possible but can lead to inconsitencies between graphviews and the parent graph's propery maps.
Here is an example of data dictionnary for the graph vertices:
vertex_data
{'height': (0.0, 'float'), 'is_active_vert': (1, 'bool'), 'is_alive': (1, 'bool'), 'is_cell_vert': (0, 'bool'), 'rho': (0.0, 'float'), 'theta': (0.0, 'float'), 'x': (0.0, 'float'), 'y': (0.0, 'float'), 'z': (0.0, 'float')}
Here we show how the hexagnal cell lattice is generated.
We go for a very explicit strategy. Note that we could use a more generic one by considering a Delaunay triangulation to a set of cell centers, followed by a Voronoï tessalation giving the junctions meshes.
print(base_grid.__doc__)
Creates a 2D hexagonal grid with the following geometry: j j j j c c c c j j j j j j j j c c c c j j j j j j j Parameters ---------- n_cells_x : int Number of cells along the first dimention (i.e the number of cell vertices on an axis parallel to the first dimention) n_cells_y : int Number of cells along the second dimention (i.e the number of cell vertices on an axis parallel to the 2nd dimention) delta_x : float, optional Spacing between the verticies along first axis, default 1 delta_y : float, optional Spacing between the verticies along 2nd axis default 1 Note ---- the total number of vertices N_v is given by (N_v = n_cells_x * n_cells_y * 3) The number of _cell_ verices is n_cells_x * n_cells_y. Returns ------- pos: np.ndarray Two dimentional array with shape `(N_v, 2)` with the vertex positions is_cell_vert: np.ndarray Array with shape `(N_v,)` with the vertex labels (1 if it is a cell vertex, 0 if it is not)
n_cells_x = 4
n_cells_y = 4
pos, is_cell_vert = base_grid(n_cells_x, n_cells_y)
fig, ax = plt.subplots(figsize=(12, 12))
ax.scatter(pos[:, 0][True - is_cell_vert],
pos[:, 1][True - is_cell_vert],
c='g', s=100)
ax.scatter(pos[:, 0][is_cell_vert],
pos[:, 1][is_cell_vert],
c='b', s=100)
ax.set_aspect('equal')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
<matplotlib.text.Text at 0x7f3627b81a58>
/home/guillaume/anaconda/envs/python3/lib/python3.4/site-packages/matplotlib-1.4.x-py3.4-linux-x86_64.egg/matplotlib/figure.py:1644: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect. warnings.warn("This figure includes Axes that are not "
Here we go from the base grid to the complete graph, by calling graph_tool geometric_graph
over the 3D positions. The base grid is projected from a plane to the
cylinder surface.
### Typical number of cells along z
n_length = 5
### Typical number of cells around the cylinder
n_circum = 7
### Distance between 2 vertices
l_0 = 1.
### heights
h_0 = 1.
### Compute the cylinder radius from the number of vertices
### (This is computed inside `cylinidrical`, but we put it there for reference)
rho_c = n_circum * l_0 / (2 * np.pi)
### Compute the lumen radius from the prefered height
rho_lumen = rho_c - h_0
delta_theta = 2 * np.pi / n_circum
delta_z = delta_theta * rho_c * np.sqrt(3)/2.
graph = cylindrical(n_length, n_circum, l_0=1, h_0=1)
2015-05-28 11:19:26,918 -leg_joint.epithelium.generation -reorient_edges -- INFO -filpped 90 edges
def draw(graph, coords, **kwargs):
'''
Draws the graph with different colors for cells and junctions
Note that the vertex and edge properties map `is_cell_vert`
and `is_junction_edge` are assumed to be present in the graph's
internal properties.
Parameters
----------
graph: :class:`gt.Graph`
coords: sequence of two strings
the names of the internal vertex property maps
used as coordinates to plot the graph
'''
is_cell_vert = graph.vp['is_cell_vert']
is_junction_edge = graph.ep['is_junction_edge']
pos_vp = gt.group_vector_property([graph.vp[c] for c in coords])
not_junction = is_junction_edge.copy()
not_junction.a = 1 - is_junction_edge.a
gt.graph_draw(graph, pos_vp, inline=True,
vertex_fill_color=is_cell_vert,
edge_color=not_junction,
output_size=(600, 400), **kwargs)
## mesh's z axis rotation around the image y axis
z_angle = np.pi/5
## Rotation of the mesh around the mesh's z axis
d_theta = np.pi/3
### Projected coordinates
graph.vp['psd_x'] = graph.vp['x'].copy()
graph.vp['psd_y'] = graph.vp['y'].copy()
pseudo_x = graph.vp['psd_x']
pseudo_y = graph.vp['psd_y']
graph.vp['rho'].a = np.hypot(graph.vp['x'].a,
graph.vp['y'].a)
pseudo_x.a = graph.vp['z'].a * np.cos(z_angle) - graph.vp['rho'].a * np.sin(
graph.vp['theta'].a + d_theta) * np.sin(z_angle)
pseudo_y.a = graph.vp['rho'].a * np.cos(graph.vp['theta'].a + d_theta)
draw(graph, ['psd_x', 'psd_y'], vorder=graph.vp['y'])
graph.vp['sigma'] = graph.vp['x'].copy()
graph.vp['sigma'].a = np.arctan2(graph.vp['y'].a, graph.vp['x'].a) * graph.vp['rho'].a
draw(graph, ['sigma', 'z'])
In order to access all the geometrical properties of the simulated tissue, we look at the graph as a mesh structure.
The unit element of this mesh is a face formed by one junction edge and two cell to junction edges, like so:
jv_i --> jv_j
^ ^
\ /
cell
By definition of the cell to junction edges, always from the cell to the junction vertex, this triangle topology is sufficient to contruct the faces.
from leg_joint.geometry.geometry import Mesh
%pdb
Automatic pdb calling has been turned ON
mesh = Mesh(graph, rho_lumen)
mesh.update_geometry()
print(sorted(mesh.cell_graph.vp.keys()))
2015-05-28 11:19:32,291 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:19:32,294 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:19:32,331 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:19:32,332 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:19:32,334 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead
['area', 'contractility', 'height', 'is_active_vert', 'is_alive', 'is_cell_vert', 'num_sides', 'perimeter', 'prefered_vol', 'psd_x', 'psd_y', 'rho', 'sigma', 'theta', 'vol', 'vol_elasticity', 'x', 'y', 'z']
mesh.update_num_sides()
mesh.cell_graph.vp['num_sides'].fa
PropertyArray([3, 3, 3, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3], dtype=int32)
mesh.graph.set_edge_filter(mesh.graph.ep['is_junction_edge'], inverted=True)
print(mesh.graph.degree_property_map('out').fa)
mesh.graph.clear_filters()
[4 0 0 4 0 0 4 0 0 4 0 0 4 0 0 4 0 0 4 0 0 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 4 0 0 4 0 0 4 0 0 4 0 0 4 0 0 4 0 0 4 0 0]
for cell in mesh.cell_graph.vertices():
if mesh.cell_graph.vp['num_sides'][cell] != mesh.graph.degree_property_map('out')[cell]:
mesh.cell_graph.vp['is_alive'][cell] = 0
print('dead cell {}'.format(cell))
dead cell 0 dead cell 3 dead cell 6 dead cell 9 dead cell 12 dead cell 15 dead cell 18 dead cell 84 dead cell 87 dead cell 90 dead cell 93 dead cell 96 dead cell 99 dead cell 102
from leg_joint.utils import _to_3d
je_idx = pd.MultiIndex.from_tuples([(mesh.junction_graph.vertex_index[s],
mesh.junction_graph.vertex_index[t])
for s, t in mesh.junction_graph.edges()],
names=('jv_i', 'jv_j'))
jv_idx = pd.Index(mesh.junction_graph.vertex_index.copy().fa, name='jv')
jv_idx
Int64Index([ 1, 2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20, 21, 22, 24, 25, 27, 28, 30, 31, 33, 34, 36, 37, 39, 40, 43, 44, 46, 47, 49, 50, 52, 53, 55, 56, 58, 59, 61, 62, 63, 64, 66, 67, 69, 70, 72, 73, 75, 76, 78, 79, 81, 82, 85, 86, 88, 89, 91, 92, 94, 95, 97, 98, 100, 101, 103, 104], dtype='int64', name='jv')
draw(mesh.junction_graph, ['sigma', 'z'])
mesh.graph.ep['dx']
<PropertyMap object with key type 'Edge' and value type 'double', for Graph 0x7f36278afe10, at 0x7f361ffe4b70, with values: [ 0.36598048 0.07355167 0.47734099 -0.03760507 0.11088188 0.31122461 -0.49372483 -0.9768541 -0.45016603 -0.94412052 0.24976694 -0.22010529 0.4939545 -0.31165804 -0.21956482 -0.49951311 0.49397527 0.1822346 -0.56122122 -0.11142236 0.98792977 -0.24999661 0.18250022 0.47723944 0.14848696 -0.45016603 0.8899737 0.07355167 -0.03760507 0.11088188 -0.22010529 0.47723944 -0.4939545 0.92740546 0.81610974 0.36618425 -0.41273426 -0.36618425 -0.9768541 0.18250022 -0.49372483 -0.41287324 -0.44992549 -0.86279873 0.41273426 0.61618085 -0.56122122 0.49951311 0.36598048 -0.7243923 0.67763852 -0.24999661 0.18250022 -0.11142236 0.61618085 -0.31122461 -0.47734099 0.24999661 0.9768541 0.86279873 0.45016603 -0.3611893 -0.41287324 -0.31165804 0.24976694 0.44992549 -0.89021423 0.36598048 0.81610974 0.61574742 0.81610974 0.41287324 0.86279873 -0.7243923 -0.36618425 -0.49397527 -0.36618425 -0.98792977 0.49951311 0.49951311 -0.99348838 0.94412052 0.47734099 0.1822346 -0.11142236 0.1822346 -0.24976694 0.98792977 0.29365696 0.47734099 -0.9768541 -0.45016603 -0.41287324 0.7243923 0.44992549 0.07355167 -0.89021423 0.4939545 -0.92740546 -0.4939545 -0.22010529 -0.49397527 0.92740546 0.31165804 0.21956482 -0.49951311 0.03733021 0.49397527 -0.99348838 0.98792977 0.94412052 -0.56122122 0.47734099 -0.47723944 -0.47723944 0.45016603 -0.8899737 0.11088188 -0.41273426 0.11142236 -0.14848696 0.47723944 -0.4939545 0.92740546 -0.4939545 -0.03733021 -0.49397527 -0.67763852 0.49372483 0.31165804 0.61574742 0.3611893 0.94412052 -0.03760507 0.41273426 -0.47723944 0.45016603 0.29365696 -0.8899737 -0.18250022 0.44992549 -0.36598048 0.41273426 -0.61574742 -0.7243923 -0.89021423 -0.11142236 0.31165804 -0.24976694 -0.61618085 -0.36598048 -0.67763852 -0.11088188 0.3611893 -0.24999661 -0.99348838 0.31122461 0.41273426 -0.31165804 -0.07355167 0.24976694 -0.29365696 0.36598048 0.61574742 -0.86279873 -0.1822346 0.11142236 -0.21956482 0.36618425 -0.24976694 0.3611893 -0.36598048 0.9768541 -0.11088188 -0.1822346 0.1822346 -0.03733021 -0.24976694 -0.31122461 0.03760507 -0.41287324 0.14848696 0.29365696 -0.1822346 -0.8899737 0.03733021 0.11142236 -0.21956482 -0.47734099 0.45016603 0.03760507 0.22010529 -0.67763852 -0.98792977 0.21956482 0.03760507 -0.11088188 0.03733021 0.18250022 -0.49372483 -0.44992549 0.14848696 -0.49397527 0.03733021 0.49951311 -0.03733021 -0.81610974 0.11088188 -0.18250022 0.56122122 -0.03760507 0.61618085 -0.14848696 -0.31122461 0.24999661 0.36618425 0.24999661 -0.49951311 -0.18250022 0.49372483 0.44992549 0.41287324 0.86279873 -0.36618425 0.41287324 -0.44992549 0.89021423 -0.24999661 0.31122461 -0.61618085 0.49397527 -0.47734099 0.99348838 -0.92740546 0.7243923 0.4939545 -0.47723944 -0.41273426 0.11142236 -0.61574742 -0.36598048 0.31165804 -0.29365696 -0.07355167 0.24976694 -0.1822346 -0.03733021 -0.31122461 0.22010529 0.03760507 -0.11088188 0.56122122 -0.81610974 -0.18250022 0.24999661 0.36618425 -0.49951311 0.89021423 0.41287324 -0.44992549 0.99348838 -0.94412052 -0.47734099 0.49397527 -0.31165804 0.4939545 0.8899737 0.47723944 -0.45016603 0.67763852 -0.41273426 -0.3611893 0.31122461]>
mesh.update_geometry()
grad_i_lij = - np.array(
[mesh.junction_graph.ep[dc].fa for dc in mesh.dcoords]).T/_to_3d(mesh.junction_graph.ep['edge_length'].fa)
grad_i_lij = pd.DataFrame(grad_i_lij, index=je_idx,
columns=mesh.coords)
grad_i_lij.head()
2015-05-28 11:41:26,398 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:41:26,401 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:41:26,440 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:41:26,442 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead 2015-05-28 11:41:26,444 -leg_joint.geometry.geometry -__getitem__ -- WARNING -Use cast instead
x | y | z | ||
---|---|---|---|---|
jv_i | jv_j | |||
1 | 2 | -0 | -0 | -0 |
22 | -0 | -0 | -0 | |
2 | 24 | -0 | -0 | -0 |
4 | 5 | -0 | -0 | -0 |
25 | -0 | -0 | -0 |
def tension_grad(mesh):
grad_t = mesh.grad_array.copy()
tensions = mesh.junction_graph.ep['line_tension'].fa
#tensions.index.names = ('jv_i', 'jv_j')
_grad_t = mesh.grad_i_lij * _to_3d(tensions)
grad_t.loc[mesh.uix_active_i] = _grad_t.sum(
level='jv_i').loc[mesh.uix_active_i].values
grad_t.loc[mesh.uix_active_j] -= _grad_t.sum(
level='jv_j').loc[mesh.uix_active_j].values
return grad_t
mesh.junction_graph.ep['line_tension'].fa = 100
tensions = mesh.junction_graph.ep['line_tension'].fa
for (e, idx) in zip(mesh.junction_graph.edges(), je_idx):
srce, trgt = idx
assert e == mesh.junction_graph.edge(srce, trgt)
_grad_t = grad_i_lij * _to_3d(tensions)
_grad_t = pd.DataFrame(_grad_t, index=je_idx)
_grad_t.head()
x | y | z | ||
---|---|---|---|---|
jv_i | jv_j | |||
1 | 2 | 56.332006 | -82.623877 | -0.000000 |
22 | -18.254275 | 46.511134 | -86.622722 | |
2 | 24 | 36.626958 | -33.984850 | -86.622722 |
4 | 5 | 99.720380 | -7.473009 | -0.000000 |
25 | -47.745223 | 14.727451 | -86.622722 |
grad_array = pd.DataFrame(np.zeros((mesh.junction_graph.num_vertices(), 3)),
index=jv_idx, columns=mesh.coords)
grad_array.head()
x | y | z | |
---|---|---|---|
jv | |||
1 | 0 | 0 | 0 |
2 | 0 | 0 | 0 |
4 | 0 | 0 | 0 |
5 | 0 | 0 | 0 |
7 | 0 | 0 | 0 |
sum_i = _grad_t.sum(level='jv_i')
sum_j = _grad_t.sum(level='jv_j')
sum_i.index.name = 'jv'
sum_j.index.name = 'jv'
grad_array.loc[sum_i.index] = sum_i
grad_array.loc[sum_j.index] -= sum_j
grad_array.head(10)
x | y | z | |
---|---|---|---|
jv | |||
1 | 38.077731 | -36.112744 | -8.662272e+01 |
2 | -19.705048 | 48.639027 | -8.662272e+01 |
4 | 51.975157 | 7.254441 | -8.662272e+01 |
5 | -50.313419 | 14.919910 | -8.662272e+01 |
7 | 26.734230 | 45.158884 | -8.662272e+01 |
8 | -43.034759 | -30.034203 | -8.662272e+01 |
10 | -18.638118 | 49.057766 | -8.662272e+01 |
11 | -14.520384 | 96.336450 | 1.421085e-14 |
13 | -49.975583 | 16.015150 | -8.662272e+01 |
14 | -73.253916 | 67.969701 | 4.263256e-14 |
cell_idx = pd.Index(mesh.cell_graph.vertex_index.copy().fa, name='cell')
a = np.arange(4)
isinstance(a, list)
False
from leg_joint.geometry.geometry import VertexFacesView, EdgeFacesView
mesh.cell_graph.vp['contractility'].fa = 100
gamma_L = mesh.cell_graph.vp['contractility'].copy()
cell_fv = VertexFacesView(mesh.cell_graph, mesh.tix_a)
gamma_L.fa = mesh.cell_graph.vp['contractility'].fa * mesh.cell_graph.vp['perimeter'].fa
mesh.cell_graph.vp['gamma_L'] = gamma_L
gamma_L = cell_fv['gamma_L']
gamma_L.index = mesh.tix_aij
area_term = gamma_L.groupby(level=('jv_i', 'jv_j')).sum()
area_term.head()
jv_i jv_j 1 2 599.161641 22 898.742461 2 24 898.742461 4 5 599.161641 25 898.742461 Name: gamma_L, dtype: float64
area_term = area_term.loc[je_idx]
_grad_c = grad_i_lij* _to_3d(area_term)
grad_c = grad_array.copy() * 0.
sum_i = _grad_c.sum(level='jv_i')
sum_j = _grad_c.sum(level='jv_j')
sum_i.index.name = 'jv'
sum_j.index.name = 'jv'
grad_c.loc[sum_i.index] = sum_i
grad_c.loc[sum_j.index] -= sum_j
for ix0, ix1 in zip(grad_i_lij.index.values,
area_term.index.values):
assert ix0 == ix1
grad_v = grad_array.copy()
grad_v[:] = 0.
mesh.cell_graph.vp['vol_elasticity'].fa = 1.
mesh.cell_graph.vp['prefered_vol'] = mesh.cell_graph.new_vertex_property('double')
mesh.cell_graph.vp['prefered_vol'].fa = 1.
mesh.cell_graph.vp['KV_V0'] = mesh.cell_graph.new_vertex_property('double')
KV_V0 = mesh.cell_graph.vp['vol_elasticity'].fa * (
mesh.cell_graph.vp['vol'].fa - mesh.cell_graph.vp['prefered_vol'].fa)
KV_V0
array([ 1.27595763, 1.27595763, 1.27595763, 1.27595763, 1.27595763, 1.27595763, 1.27595763, 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 7.1399895 , 1.27595763, 1.27595763, 1.27595763, 1.27595763, 1.27595763, 1.27595763, 1.27595763])
def volume_grad(mesh):
'''
Computes :math:`\sum_\alpha\nabla_i \left(K (V_\alpha - V_0)^2\right)`
'''
grad_v = mesh.grad_array.copy()
grad_v[:] = 0
elasticity = mesh.udf_cell['vol_elasticity']
pref_V = mesh.udf_cell['prefered_vol']
V = mesh.udf_cell['vol']
KV_V0 = elasticity * (V - pref_V)
tri_KV_V0 = KV_V0.loc[mesh.tix_a]
tri_KV_V0.index = mesh.tix_aij
r_ijs = mesh.tdf_itoj[mesh.dcoords]
cross_ur = pd.DataFrame(np.cross(mesh.faces[mesh.normal_coords], r_ijs),
index=mesh.tix_aij, columns=mesh.coords)
h_nu = mesh.udf_cell['height'] / (2 * mesh.udf_cell['num_sides'])
grad_i_V_cell = cross_ur.sum(level='cell') * _to_3d(KV_V0 * h_nu)
cell_term_i = grad_i_V_cell.loc[mesh.tix_a].set_index(mesh.tix_ai)
cell_term_j = grad_i_V_cell.loc[mesh.tix_a].set_index(mesh.tix_aj)
grad_v.loc[mesh.uix_active_i] += cell_term_i.loc[mesh.uix_ai].sum(
level='jv_i').loc[mesh.uix_active_i].values/2
grad_v.loc[mesh.uix_active_j] += cell_term_j.loc[mesh.uix_aj].sum(
level='jv_j').loc[mesh.uix_active_j].values/2
_r_to_rho_i = mesh.udf_jv_i[mesh.coords] / _to_3d(mesh.udf_jv_i['rho'])
_r_to_rho_j = mesh.udf_jv_j[mesh.coords] / _to_3d(mesh.udf_jv_j['rho'])
r_to_rho_i = _r_to_rho_i.loc[mesh.tix_i].set_index(mesh.tix_aij)
r_to_rho_j = _r_to_rho_j.loc[mesh.tix_j].set_index(mesh.tix_aij)
r_ai = mesh.tdf_atoi[mesh.dcoords]
r_aj = mesh.tdf_atoj[mesh.dcoords]
normals = mesh.faces[mesh.normal_coords]
cross_ai = pd.DataFrame(np.cross(normals, r_ai),
index=mesh.tix_aij, columns=mesh.coords)
cross_aj = pd.DataFrame(np.cross(normals, r_aj),
index=mesh.tix_aij, columns=mesh.coords)
tri_height = mesh.tdf_cell['height']
tri_height.index = mesh.tix_aij
sub_area = mesh.faces['sub_area']
_ij_term = _to_3d(tri_KV_V0) *(_to_3d(sub_area / 2) * r_to_rho_i
- _to_3d(tri_height / 2) * cross_aj)
_jk_term = _to_3d(tri_KV_V0) *(_to_3d(sub_area / 2) * r_to_rho_j
+ _to_3d(tri_height / 2) * cross_ai)
#ij_term = _ij_term.groupby(level=('jv_i', 'jv_j')).sum()
#jk_term = _jk_term.groupby(level=('jv_j', 'jv_i')).sum()
grad_v.loc[mesh.uix_active_i] += _ij_term.sum(
level='jv_i').loc[mesh.uix_active_i].values
grad_v.loc[mesh.uix_active_j] += _jk_term.sum(
level='jv_j').loc[mesh.uix_active_j].values
return grad_v