volume of a irregular tetrahedron
426463
def vol_tehrahedron(poly):
    """"""
    p_a = np.array(poly[0])
    p_b = np.array(poly[1])
    p_c = np.array(poly[2])
    p_d = np.array(poly[3])
    return abs(np.dot(
        np.subtract(p_a, p_d),
        np.cross(
            np.subtract(p_b, p_d),
            np.subtract(p_c, p_d))) / 6)
1346947
def tricu(P, k=0):
    """"""
    tri = numpy.sum(numpy.mgrid[[slice(0,_,1) for _ in P.shape]], 0)
    tri = tri<len(tri) + k

    if isinstance(P, Poly):
        A = P.A.copy()
        B = {}
        for key in P.keys:
            B[key] = A[key]*tri
        return Poly(B, shape=P.shape, dim=P.dim, dtype=P.dtype)

    out = P*tri
    return out
575058
def irregular_contour(x, y, z, func=plt.contourf, func_kwargs=dict(),
                      grid_size=(100,100), padding_fraction=0.05,
                      interp_method='nearest'):
  '''
  '''
  from scipy.interpolate import griddata  # Late import; scipy is optional
  x, y, z = map(np.asanyarray, (x, y, z))
  x_range = (x.min(), x.max())
  y_range = (y.min(), y.max())
  pad_x = padding_fraction * -np.subtract.reduce(x_range)
  pad_y = padding_fraction * -np.subtract.reduce(y_range)
  grid_x = np.linspace(x_range[0] - pad_x, x_range[1] + pad_x, grid_size[0])
  grid_y = np.linspace(y_range[0] - pad_y, y_range[1] + pad_y, grid_size[1])
  grid_z = griddata((x, y), z, (grid_x[None], grid_y[:,None]),
                    method=interp_method)
  return func(grid_x, grid_y, grid_z, **func_kwargs)
520439
def euclidean3d(v1, v2):
    """"""
    if not len(v1) == 3 and len(v2) == 3:
        print("Vectors are not in 3D space. Returning None.")
        return None
    return np.sqrt((v1[0] - v2[0]) ** 2 + (v1[1] - v2[1]) ** 2 + (v1[2] - v2[2]) ** 2)
244464
def near_reduce(self, coords_set, threshold=1e-4):
        
        unique_coords = []
        coords_set = [self.slab.lattice.get_fractional_coords(coords)
                      for coords in coords_set]
        for coord in coords_set:
            if not in_coord_list_pbc(unique_coords, coord, threshold):
                unique_coords += [coord]
        return [self.slab.lattice.get_cartesian_coords(coords)
                for coords in unique_coords]
471561
def partsphere(self, x):
        """"""
        self.counter += 1
        # return np.random.rand(1)[0]**0 * sum(x**2) + 1 * np.random.rand(1)[0]
        dim = len(x)
        x = array([x[i % dim] for i in xrange(2 * dim)])
        N = 8
        i = self.counter % dim
        # f = sum(x[i:i + N]**2)
        f = sum(x[np.random.randint(dim, size=N)]**2)
        return f
1285
def triangulate_quads(quads):
    
    if len(quads) == 0:
        return quads
    quads = np.asanyarray(quads)
    faces = np.vstack((quads[:, [0, 1, 2]],
                       quads[:, [2, 3, 0]]))
    return faces
1117930
def fill_volume(volume)
      
      # First fill a cube that bounds the volume
      max = volume.max_point
      min = volume.min_point
      
      dx = max[0] - min[0]
      dy = max[1] - min[1]
      dz = max[2] - min[2]
      
      bulk = get_bulk
      
      # This inverse matrix gives the number of repetitions 
      m = Matrix[[dx,0,0], [0,dy,0], [0,0,dz]]
      v = Matrix[bulk.lattice_vectors[0].to_a, 
                 bulk.lattice_vectors[1].to_a,
                 bulk.lattice_vectors[2].to_a]
      rep_mat = m*(v.inverse)
      
      # The only way I can figure out how to do this for an 
      # arbitrary set of lattice vectors is to fill the volume
      # out along each edge of the super-cube and then eliminate duplicates
      atoms = []
      3.times do |i|
        # this vector is the number of repetitions in the unit cell
        # to fill the volume out along the i-th edge of the super-cube
        n_repeat = rep_mat.row(i)
        
        # Give the proper sign to the repeat
        nx = (n_repeat[0] < 0) ? n_repeat[0].floor-1 : n_repeat[0].ceil+1
        ny = (n_repeat[1] < 0) ? n_repeat[1].floor-1 : n_repeat[1].ceil+1
        nz = (n_repeat[2] < 0) ? n_repeat[2].floor-1 : n_repeat[2].ceil+1
        
        atoms += bulk.repeat(nx, ny, nz).atoms.find_all{|a| volume.contains_point(a.x, a.y, a.z)}
      end
      Geometry.new(atoms.uniq)
    end
987795
def _compute_cell_circumcenters(self):
        
        # Just like for triangular cells, tetrahedron circumcenters are most easily
        # computed with the quadrilateral coordinates available.
        # Luckily, we have the circumcenter-face distances (cfd):
        #
        #   CC = (
        #       + cfd[0] * face_area[0] / sum(cfd*face_area) * X[0]
        #       + cfd[1] * face_area[1] / sum(cfd*face_area) * X[1]
        #       + cfd[2] * face_area[2] / sum(cfd*face_area) * X[2]
        #       + cfd[3] * face_area[3] / sum(cfd*face_area) * X[3]
        #       )
        #
        # (Compare with
        # <https://en.wikipedia.org/wiki/Trilinear_coordinates#Between_Cartesian_and_trilinear_coordinates>.)
        # Because of
        #
        #    cfd = zeta / (24.0 * face_areas) / self.cell_volumes[None]
        #
        # we have
        #
        #   CC = sum_k (zeta[k] / sum(zeta) * X[k]).
        #
        alpha = self._zeta / numpy.sum(self._zeta, axis=0)

        self._circumcenters = numpy.sum(
            alpha[None].T * self.node_coords[self.cells["nodes"]], axis=1
        )
        return
664268
def lattice(self, lattice):
        
        self._lattice = lattice
        self._coords = self._lattice.get_cartesian_coords(self._frac_coords)
507002
def volume( self ):
        
        return np.dot( self.matrix[0], np.cross( self.matrix[1], self.matrix[2] ) )
608638
def _g_3(self):
        """"""
        # return 3 * (1.0 - self._n_3()) / (self._vertices_omegas[3] - self._omega)
        return (3 * self._f(1, 3) * self._f(2, 3) /
                (self._vertices_omegas[3] - self._vertices_omegas[0]))
367278
def get_grid_dimensions(current_problem_size, params, grid_div, block_size_names):
    """"""
    def get_dimension_divisor(divisor_list, default, params):
        if divisor_list is None:
            if default in params:
                divisor_list = [default]
            else:
                return 1
        return numpy.prod([int(eval(replace_param_occurrences(s, params))) for s in divisor_list])
    divisors = [get_dimension_divisor(d, block_size_names[i], params) for i, d in enumerate(grid_div)]
    return tuple(int(numpy.ceil(float(current_problem_size[i]) / float(d))) for i, d in enumerate(divisors))
959573
def mass_3d(self, r, rho0, Ra, Rs):
        
        m_3d = 4 * np.pi * rho0 * Ra ** 2 * Rs ** 2 / (Rs ** 2 - Ra ** 2) * (Rs * np.arctan(r / Rs) - Ra * np.arctan(r / Ra))
        return m_3d
147824
function ( radius, detail ) {

	var vertices = [
		1, 1, 1, 	- 1, - 1, 1, 	- 1, 1, - 1, 	1, - 1, - 1
	];

	var indices = [
		2, 1, 0, 	0, 3, 2,	1, 3, 0,	2, 3, 1
	];

	PolyhedronBufferGeometry.call( this, vertices, indices, radius, detail );

	this.type = '';

	this.parameters = {
		radius: radius,
		detail: detail
	};

}
115116
def noise3(self, x, y, z):
		
		# Skew the input space to determine which simplex cell we're in
		s = (x + y + z) * _F3
		i = floor(x + s)
		j = floor(y + s)
		k = floor(z + s)
		t = (i + j + k) * _G3
		x0 = x - (i - t) # "Unskewed" distances from cell origin
		y0 = y - (j - t)
		z0 = z - (k - t)

		# For the 3D case, the simplex shape is a slightly irregular tetrahedron. 
		# Determine which simplex we are in. 
		if x0 >= y0:
			if y0 >= z0:
				i1 = 1; j1 = 0; k1 = 0
				i2 = 1; j2 = 1; k2 = 0
			elif x0 >= z0:
				i1 = 1; j1 = 0; k1 = 0
				i2 = 1; j2 = 0; k2 = 1
			else:
				i1 = 0; j1 = 0; k1 = 1
				i2 = 1; j2 = 0; k2 = 1
		else: # x0 < y0
			if y0 < z0:
				i1 = 0; j1 = 0; k1 = 1
				i2 = 0; j2 = 1; k2 = 1
			elif x0 < z0:
				i1 = 0; j1 = 1; k1 = 0
				i2 = 0; j2 = 1; k2 = 1
			else:
				i1 = 0; j1 = 1; k1 = 0
				i2 = 1; j2 = 1; k2 = 0
		
		# Offsets for remaining corners
		x1 = x0 - i1 + _G3
		y1 = y0 - j1 + _G3
		z1 = z0 - k1 + _G3
		x2 = x0 - i2 + 2.0 * _G3
		y2 = y0 - j2 + 2.0 * _G3
		z2 = z0 - k2 + 2.0 * _G3
		x3 = x0 - 1.0 + 3.0 * _G3
		y3 = y0 - 1.0 + 3.0 * _G3
		z3 = z0 - 1.0 + 3.0 * _G3

		# Calculate the hashed gradient indices of the four simplex corners
		perm = self.permutation
		ii = int(i) % self.period
		jj = int(j) % self.period
		kk = int(k) % self.period
		gi0 = perm[ii + perm[jj + perm[kk]]] % 12
		gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12
		gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12
		gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12

		# Calculate the contribution from the four corners
		noise = 0.0
		tt = 0.6 - x0**2 - y0**2 - z0**2
		if tt > 0:
			g = _GRAD3[gi0]
			noise = tt**4 * (g[0] * x0 + g[1] * y0 + g[2] * z0)
		else:
			noise = 0.0
		
		tt = 0.6 - x1**2 - y1**2 - z1**2
		if tt > 0:
			g = _GRAD3[gi1]
			noise += tt**4 * (g[0] * x1 + g[1] * y1 + g[2] * z1)
		
		tt = 0.6 - x2**2 - y2**2 - z2**2
		if tt > 0:
			g = _GRAD3[gi2]
			noise += tt**4 * (g[0] * x2 + g[1] * y2 + g[2] * z2)
		
		tt = 0.6 - x3**2 - y3**2 - z3**2
		if tt > 0:
			g = _GRAD3[gi3]
			noise += tt**4 * (g[0] * x3 + g[1] * y3 + g[2] * z3)
		
		return noise * 32.0
291673
private void tred2()
    {
        for(int j = 0; j < n; j++)
            d[j] = V.get(n-1, j);

        // Householder reduction to tridiagonal form.

        for (int i = n - 1; i > 0; i--)
        {

            // Scale to avoid under/overflow.

            double scale = 0.0;
            double h = 0.0;
            for (int k = 0; k < i; k++)
            {
                scale = scale + abs(d[k]);
            }
            if (scale == 0.0)
            {
                e[i] = d[i - 1];
                for (int j = 0; j < i; j++)
                {
                    d[j] = V.get(i-1, j);
                    V.set(i, j, 0.0);
                    V.set(j, i, 0.0);
                }
            }
            else
            {

                // Generate Householder vector.

                for (int k = 0; k < i; k++)
                {
                    d[k] /= scale;
                    h += d[k] * d[k];
                }
                double f = d[i - 1];
                double g = sqrt(h);
                if (f > 0)
                    g = -g;
                e[i] = scale * g;
                h -= f * g;
                d[i - 1] = f - g;
                Arrays.fill(e, 0, i, 0.0);

                // Apply similarity transformation to remaining columns.

                for (int j = 0; j < i; j++)
                {
                    f = d[j];
                    V.set(j, i, f);
                    g = e[j] + V.get(j, j) * f;
                    for (int k = j + 1; k <= i - 1; k++)
                    {
                        g += V.get(k,j) * d[k];
                        e[k] += V.get(k, j) * f;
                    }
                    e[j] = g;
                }
                f = 0.0;
                for (int j = 0; j < i; j++)
                {
                    e[j] /= h;
                    f += e[j] * d[j];
                }
                double hh = f / (h + h);
                for (int j = 0; j < i; j++)
                {
                    e[j] -= hh * d[j];
                }
                for (int j = 0; j < i; j++)
                {
                    f = d[j];
                    g = e[j];
                    
                    for (int k = j; k <= i - 1; k++)
                    {
                        V.increment(k, j, -(f * e[k] + g * d[k]));
                    }
                    d[j] = V.get(i-1, j);
                    V.set(i, j, 0.0);
                }
            }
            d[i] = h;
        }


        // Accumulate transformations.

        for (int i = 0; i < n - 1; i++)
        {
            V.set(n-1, i, V.get(i, i));
            V.set(i, i, 1.0);
            double h = d[i + 1];
            if (h != 0.0)
            {
                for (int k = 0; k <= i; k++)
                {
                    d[k] = V.get(k, i+1) / h;
                }
                for (int j = 0; j <= i; j++)
                {
                    double g = 0.0;
                    for (int k = 0; k <= i; k++)
                    {
                        g += V.get(k, i+1) * V.get(k, j);
                    }
                    
                    RowColumnOps.addMultCol(V, j, 0, i+1, -g, d);
                }
            }

            RowColumnOps.fillCol(V, i+1, 0, i+1, 0.0);
            
        }
        for (int j = 0; j < n; j++)
        {
            d[j] = V.get(n-1, j);
            V.set(n-1, j, 0.0);
        }
        V.set(n-1, n-1, 1.0);
        e[0] = 0.0;
    }