Commit b7378502 authored by Germain Clavier's avatar Germain Clavier
Browse files

Several modifications to the following files to make the new examples work.

Plotview is annoying. Duck it. Duck all visualisation libs form this repo.
parent 4809ebef
Loading
Loading
Loading
Loading
+3 −1
Original line number Diff line number Diff line
# simple test of animate tool
# requires files/bucky*png files

import animate

a = animate("files/bucky*gif")
a.play()
a.delay(0.1)
@@ -11,4 +13,4 @@ a.next()
a.previous()
a.frame(1)

print "all done ... type CTRL-D to exit Pizza.py"
print("all done ... type CTRL-D to exit Pizza.py")
+13 −11
Original line number Diff line number Diff line
@@ -67,9 +67,9 @@ class animate:

        # convert filestr into full list of files

        list = str.split(filestr)
        files = str.split(filestr)
        self.files = []
        for file in list:
        for file in files:
            self.files += glob.glob(file)
        self.nframes = len(self.files)
        if self.nframes == 0:
@@ -79,22 +79,20 @@ class animate:

        # load all images

        tkroot = Tk()
        self.tkroot = tkroot

        win1 = Toplevel(tkroot)
        win1.title("Pizza.py animate tool")

        self.images = []
        for i in range(self.nframes):
            self.images.append(PhotoImage(file=self.files[i]))

        # grab Tk instance from main

        # from __main__ import tkroot

        tkroot = Tk()
        self.tkroot = tkroot

        # GUI control window

        win1 = Toplevel(tkroot)
        win1.title("Pizza.py animate tool")

        holder1 = Frame(win1)
        button1 = Button(holder1, text="<<", command=self.first).pack(side=LEFT)
        button2 = Button(holder1, text="<", command=self.previous).pack(side=LEFT)
@@ -137,13 +135,17 @@ class animate:
        win2 = Toplevel(tkroot)
        self.image_pane = Label(win2, image=self.images[0])
        self.image_pane.pack(side=BOTTOM)
        tkroot.update_idletasks()  # force window to appear
        self.tkroot.update_idletasks()  # force window to appear

        # display 1st image

        self.index = 0
        self.display(self.index)

        # More comfortable to have a small delay IMO
        self.delay(0.1)
        self.tkroot.mainloop()

    # --------------------------------------------------------------------

    def first(self):
+50 −24
Original line number Diff line number Diff line
@@ -160,6 +160,7 @@ import sys
import glob
from os import popen
from math import sqrt, pi, cos, sin, fabs
import numpy as np
from copy import deepcopy

try:
@@ -547,29 +548,39 @@ class cdata:

        cmd1 = teststr.replace("$x", "o.vertices[v1][0]")
        cmd1 = cmd1.replace("$y", "o.vertices[v1][1]")
        cmd1 = "flag1 = " + cmd1.replace("$z", "o.vertices[v1][2]")
        ccmd1 = compile(cmd1, "", "single")
        # cmd1 = "flag1 = " + cmd1.replace("$z", "o.vertices[v1][2]")
        cmd1 = cmd1.replace("$z", "o.vertices[v1][2]")
        # ccmd1 = compile(cmd1, "", "single")

        cmd2 = teststr.replace("$x", "o.vertices[v2][0]")
        cmd2 = cmd2.replace("$y", "o.vertices[v2][1]")
        cmd2 = "flag2 = " + cmd2.replace("$z", "o.vertices[v2][2]")
        ccmd2 = compile(cmd2, "", "single")
        # cmd2 = "flag2 = " + cmd2.replace("$z", "o.vertices[v2][2]")
        cmd2 = cmd2.replace("$z", "o.vertices[v2][2]")
        # ccmd2 = compile(cmd2, "", "single")

        cmd3 = teststr.replace("$x", "o.vertices[v3][0]")
        cmd3 = cmd3.replace("$y", "o.vertices[v3][1]")
        cmd3 = "flag3 = " + cmd3.replace("$z", "o.vertices[v3][2]")
        ccmd3 = compile(cmd3, "", "single")
        # cmd3 = "flag3 = " + cmd3.replace("$z", "o.vertices[v3][2]")
        cmd3 = cmd3.replace("$z", "o.vertices[v3][2]")
        # ccmd3 = compile(cmd3, "", "single")
        print(cmd1)
        print(cmd2)
        print(cmd3)

        # loop over triangles in id_surf
        # 3 vertices must satisfy all 3 tests for tri's inclusion in new surf obj
        # 3 vertices must satisfy all 3 tests
        # for tri's inclusion in new surf obj

        for tri in o.triangles:
            v1 = tri[0] - 1
            v2 = tri[1] - 1
            v3 = tri[2] - 1
            exec(ccmd1)
            exec(ccmd2)
            exec(ccmd3)
            # flag1 = o.vertices[v1][2] < 2.0
            # flag2 = o.vertices[v2][2] < 2.0
            # flag3 = o.vertices[v3][2] < 2.0
            flag1 = eval(cmd1)
            flag2 = eval(cmd2)
            flag3 = eval(cmd3)
            # Flags defined in commands
            if flag1 and flag2 and flag3:
                obj.vertices.append(o.vertices[v1][:])
@@ -616,7 +627,8 @@ class cdata:
        if out_id:
            out_obj = self.objs[self.ids[out_id]]

        # pre-process SURFACE objects to bin their triangles for faster searching
        # pre-process SURFACE objects to bin
        # their triangles for faster searching

        if in_obj.style == SURFACE:
            in_obj.inside_prep()
@@ -629,15 +641,23 @@ class cdata:
        xsize = xhi - xlo
        ysize = yhi - ylo
        zsize = zhi - zlo
        print("bbox:")
        print(xlo, xhi, xsize)
        print(ylo, yhi, ysize)
        print(zlo, zhi, zsize)

        # generate particles until have enough that satisfy in/out constraints

        count = attempt = 0
        while count < npart:
            attempt += 1
            x = xlo + self.random() * xsize
            y = ylo + self.random() * ysize
            z = zlo + self.random() * zsize
            x = xlo + np.random.random() * xsize
            y = ylo + np.random.random() * ysize
            z = zlo + np.random.random() * zsize
            # x = xlo + self.random() * xsize
            # y = ylo + self.random() * ysize
            # z = zlo + self.random() * zsize
            # print(x, y, z)
            if not in_obj.inside(x, y, z):
                continue
            if out_id and out_obj.inside(x, y, z):
@@ -666,7 +686,9 @@ class cdata:
        obj.xyz = []

        on_obj = self.objs[self.ids[on_id]]
        if on_obj.style != SURFACE and on_obj.style != REGION and on_obj.style != UNION:
        if (on_obj.style != SURFACE
           and on_obj.style != REGION
           and on_obj.style != UNION):
            sys.exit("Illegal ID to place particles on")
        totalarea = on_obj.area()

@@ -1478,7 +1500,8 @@ class Surface:
    # project 3d triangles to xy plane
    # only need examine triangles in bin that point is in
    # if pt is outside bins, don't need to check
    # x,y,z is inside surf if line segment intersects an odd number of triangles
    # x,y,z is inside surf if line segment intersects
    # an odd number of triangles
    # intersection test:
    #   is xy pt in bounding rectangle of tri in xy plane ?
    #   is xy pt inside tri in xy plane (including edges and vertices) ?
@@ -1488,6 +1511,7 @@ class Surface:
    def inside(self, x, y, z):
        ix = int((x - self.xlo) * self.dxinv)
        iy = int((y - self.ylo) * self.dyinv)
        print(ix, self.nbinx, iy, self.nbiny)
        if ix < 0 or ix >= self.nbinx or iy < 0 or iy >= self.nbiny:
            return 0
        n = len(self.bin[ix][iy])
@@ -1513,7 +1537,8 @@ class Surface:
                continue

            # is x,y inside 2d triangle ?
            # cross product of each edge with vertex-to-point must have same sign
            # cross product of each edge
            # with vertex-to-point must have same sign
            # cross product = 0 is OK, means point is on edge or vertex

            c1 = (v2[0]-v1[0]) * (y-v1[1]) - (v2[1]-v1[1]) * (x-v1[0])
@@ -1535,7 +1560,8 @@ class Surface:
            # p = (x,y) - v1, v = v2 - v1, w = v3 - v2
            # find alpha,beta such that p = alpha v + beta w
            # by solving system of 2 linear eqs for their intersection
            # denom (vx - vy * wx/xy) cannot be 0 since would imply vx/vy = wx/wy
            # denom (vx - vy * wx/xy) cannot be 0
            # since would imply vx/vy = wx/wy
            # and thus e1 would be parallel to e2
            # if wy = 0, vy cannot be 0, since e1 would be parallel to e2

+29 −17
Original line number Diff line number Diff line
@@ -54,6 +54,7 @@ c.write("data.file") write out all built chains to LAMMPS data file

import math
import sys
import numpy as np
from data import data

# Class definition
@@ -74,6 +75,7 @@ class chain:
        self.seed = 12345
        self.mtype = 1
        self.btype = 1
        self.charge = 0.
        self.blen = 0.97
        self.dmin = 1.02
        self.id = "chain"
@@ -92,7 +94,9 @@ class chain:
        self.zlo = -self.zprd / 2.0
        self.zhi = self.zprd / 2.0

        print("Simulation box: %g by %g by %g" % (self.xprd, self.yprd, self.zprd))
        print("Simulation box: %g by %g by %g" %
              (self.xprd, self.yprd, self.zprd)
              )

    # --------------------------------------------------------------------

@@ -101,26 +105,29 @@ class chain:
            atoms = []
            bonds = []
            id_atom_prev = id_mol_prev = id_bond_prev = 0
            if len(self.atoms):
            if self.atoms:
                id_atom_prev = self.atoms[-1][0]
                id_mol_prev = self.atoms[-1][1]
            if len(self.bonds):
            if self.bonds:
                id_bond_prev = self.bonds[-1][0]

            for imonomer in range(nper):
                if imonomer == 0:
                    x = self.xlo + self.random() * self.xprd
                    y = self.ylo + self.random() * self.yprd
                    z = self.zlo + self.random() * self.zprd
                    # x = self.xlo + self.random() * self.xprd
                    # y = self.ylo + self.random() * self.yprd
                    # z = self.zlo + self.random() * self.zprd
                    x = np.random.uniform(self.xlo, self.xhi)
                    y = np.random.uniform(self.xlo, self.xhi)
                    z = np.random.uniform(self.xlo, self.xhi)
                    ix = iy = iz = 0
                else:
                    restriction = True
                    while restriction:
                        rsq = 2.0
                        while rsq > 1.0:
                            dx = 2.0 * self.random() - 1.0
                            dy = 2.0 * self.random() - 1.0
                            dz = 2.0 * self.random() - 1.0
                            dx = np.random.uniform(-1, 1)
                            dy = np.random.uniform(-1, 1)
                            dz = np.random.uniform(-1, 1)
                            rsq = dx * dx + dy * dy + dz * dz
                        r = math.sqrt(rsq)
                        dx, dy, dz = dx / r, dy / r, dz / r
@@ -148,7 +155,11 @@ class chain:
                else:
                    sys.exit("chain ID is not a valid value")

                atoms.append([idatom, idmol, self.mtype, x, y, z, ix, iy, iz])
                atcharge = self.charge
                atoms.append(
                        [idatom, idmol, self.mtype, atcharge,
                         x, y, z, ix, iy, iz]
                        )
                if imonomer:
                    bondid = id_bond_prev + imonomer
                    bonds.append([bondid, self.btype, idatom - 1, idatom])
@@ -188,12 +199,12 @@ class chain:

        lines = []
        for i in range(atypes):
            lines.append("%d 1.0\n" % (i + 1))
            lines.append("%d 1.0" % (i + 1))
        d.sections["Masses"] = lines

        lines = []
        for atom in self.atoms:
            line = "%d %d %d %g %g %g %d %d %d\n" % (
            line = "{:>} {:>} {:>} {:>} {:>} {:>} {:>} {:>} {:>} {:>}".format(
                atom[0],
                atom[1],
                atom[2],
@@ -203,13 +214,14 @@ class chain:
                atom[6],
                atom[7],
                atom[8],
                atom[9]
            )
            lines.append(line)
        d.sections["Atoms"] = lines

        lines = []
        for bond in self.bonds:
            line = "%d %d %d %d\n" % (bond[0], bond[1], bond[2], bond[3])
            line = "%d %d %d %d" % (bond[0], bond[1], bond[2], bond[3])
            lines.append(line)
        d.sections["Bonds"] = lines

+2 −1
Original line number Diff line number Diff line
@@ -277,11 +277,12 @@ class data:
                        )
                    )
                else:
                    f.write("{}\n".format(self.headers[keyword], keyword))
                    f.write("{} {}\n".format(self.headers[keyword], keyword))
        for pair in skeywords:
            keyword = pair[0]
            if keyword in self.sections:
                f.write("\n%s\n" % keyword)
                f.write("\n")
                for line in self.sections[keyword]:
                    f.write("{}\n".format(line))
        f.close()
Loading