Python simple.regularGrid函数代码示例

def centerline(F,dir,nx=2,mode=2,th=0.2):
    """Compute the centerline in the direction dir.

    """
    bb = F.bbox()
    x0 = F.center()
    x1 = F.center()
    x0[dir] = bb[0][dir]
    x1[dir] = bb[1][dir]
    n = array((0,0,0))
    n[dir] = nx
    
    grid = simple.regularGrid(x0,x1,n).reshape((-1,3))

    if mode > 0:
        th *= (x1[dir]-x0[dir])/nx
        n = zeros((3,))
        n[dir] = 1.0
        center = []
        for P in grid:
            test = abs(F.distanceFromPlane(P,n)) < th
            if mode == 1:
                C = F.coords[test].mean(axis=0)
            elif mode == 2:
                test = test.sum(axis=-1)
                G = F.select(test==F.coords.shape[1])
                print(G)
                C = G.center()
            center.append(C)
        grid = array(center)

    return Formex(connectPoints(grid))

def drawGridPlanes(x0,x1,nx):
    """Draw a 3D rectangular grid of planes.
        
    A grid of planes parallel to the axes is drawn in the domain bounded
    by the rectangular box [x0,x1]. The grid has nx divisions in the axis
    directions, thus planes will be drawn at nx[i]+1 positions in direction i.
    If nx[i] == 0, planes are only drawn for the initial coordinate x0.
    Thus nx=(0,2,3) results in a grid of 3x4 planes // x and
    one plane // (y,z) at coordinate x=x0[0].
    """
    x0 = asarray(x0)
    x1 = asarray(x1)
    nx = asarray(nx)

    for i in range(3):
        axes = (asarray([1,2]) + i) % 3
        if all(nx[axes] > 0):
            j,k = axes
            base = simple.regularGrid(x0[i],x1[i],nx[i]).ravel()
            x = zeros((base.shape[0],4,3))
            corners = array([x0[axes],[x1[j],x0[k]],x1[axes],[x0[j],x1[k]]])
            for j in range(4):
                x[:,j,i] = base
            x[:,:,axes] = corners
            GL.glBegin(GL.GL_QUADS)
            for p in x.reshape((-1,3)):
                GL.glVertex3fv(p)
            GL.glEnd()

def create():
    """Create a closed surface and a set of points."""
    nx,ny,nz = npts

    # Create surface
    if surface == 'file':
        S = TriSurface.read(filename).centered()
    elif surface == 'sphere':
        S = simple.sphere(ndiv=grade)

    if refine > S.nedges():
        S = S.refine(refine)

    draw(S, color='red')

    if not S.isClosedManifold():
        warning("This is not a closed manifold surface. Try another.")
        return None,None

    # Create points

    if points == 'grid':
        P = simple.regularGrid([-1.,-1.,-1.],[1., 1., 1.],[nx-1,ny-1,nz-1])
    else:
        P = random.rand(nx*ny*nz*3)

    sc = array(scale)
    siz = array(S.sizes())
    tr = array(trl)
    P = Formex(P.reshape(-1, 3)).resized(sc*siz).centered().translate(tr*siz)
    draw(P, marksize=1, color='black')
    zoomAll()

    return S,P

def drawGridLines(x0,x1,nx):
    """Draw a 3D rectangular grid of lines.
        
    A grid of lines parallel to the axes is drawn in the domain bounded
    by the rectangular box [x0,x1]. The grid has nx divisions in the axis
    directions, thus lines will be drawn at nx[i]+1 positions in direction i.
    If nx[i] == 0, lines are only drawn for the initial coordinate x0.
    Thus nx=(0,2,3) results in a grid of 3x4 lines in the plane // (y,z) at
    coordinate x=x0[0].
    """
    x0 = asarray(x0)
    x1 = asarray(x1)
    nx = asarray(nx)

    for i in range(3):
        if nx[i] > 0:
            axes = (asarray([1,2]) + i) % 3
            base = simple.regularGrid(x0[axes],x1[axes],nx[axes]).reshape((-1,2))
            x = zeros((base.shape[0],2,3))
            x[:,0,axes] = base
            x[:,1,axes] = base
            x[:,0,i] = x0[i]
            x[:,1,i] = x1[i]
            GL.glBegin(GL.GL_LINES)
            for p in x.reshape((-1,3)):
                GL.glVertex3fv(p)
            GL.glEnd()

def centerline(F,dir,nx=2,mode=2,th=0.2):
    """Compute the centerline in the direction dir.

    """
    bb = F.bbox()
    x0 = F.center()
    x1 = F.center()
    x0[dir] = bb[0][dir]
    x1[dir] = bb[1][dir]
    n = array((0,0,0))
    n[dir] = nx
    
    grid = simple.regularGrid(x0,x1,n).reshape((-1,3))

    th *= (x1[dir]-x0[dir])/nx
    n = zeros((3,))
    n[dir] = 1.0

    def localCenter(X,P,n):
        """Return the local center of points X in the plane P,n"""
        test = abs(X.distanceFromPlane(P,n)) < th # points close to plane
        if mode == 1:
            C = X[test].center()
        elif mode == 2:
            C = X[test].centroid()
        return C

    center = [ localCenter(F,P,n) for P in grid ]
    return PolyLine(center)

    def setCamera(self,bbox=None,angles=None):
        """Sets the camera looking under angles at bbox.

        This function sets the camera angles and adjusts the zooming.
        The camera distance remains unchanged.
        
        If a bbox is specified, the camera will be zoomed to make the whole
        bbox visible.
        If no bbox is specified, the current scene bbox will be used.
        If no current bbox has been set, it will be calculated as the
        bbox of the whole scene.

        If no camera angles are given, the camera orientation is kept.
        angles can be a set of 3 angles, or a string
        """
        self.makeCurrent()
        # go to a distance to have a good view with a 45 degree angle lens
        if bbox is not None:
            self.setBbox(bbox)
        #GD.debug("USING BBOX: %s" % self.bbox)
        X0,X1 = self.bbox
        center = 0.5*(X0+X1)
        # calculating the bounding circle: this is rather conservative
        self.camera.setCenter(*center)
        if type(angles) is str:
            angles = self.view_angles.get(angles)
        if angles is not None:
            try:
                self.camera.setAngles(angles)
            except:
                raise ValueError,'Invalid view angles specified'
        # Currently, we keep the default fovy/aspect
        # and change the camera distance to focus
        fovy = self.camera.fovy
        #GD.debug("FOVY: %s" % fovy)
        self.camera.setLens(fovy,self.aspect)
        # Default correction is sqrt(3)
        correction = float(GD.cfg.get('gui/autozoomfactor',1.732))
        tf = tand(fovy/2.)

        import simple,coords
        bbix = simple.regularGrid(X0,X1,[1,1,1])
        bbix = dot(bbix,self.camera.rot[:3,:3])
        bbox = coords.Coords(bbix).bbox()
        dx,dy = bbox[1][:2] - bbox[0][:2]
        hsize = max(dx,dy/self.aspect)
        offset = abs(bbox[1][2]+bbox[0][2])
        #print "hsize,offset = %s,%s" % (hsize,offset)
        dist = (hsize/tf + offset) / correction
        #print "new dist = %s" % (dist)
        
        if dist == nan or dist == inf:
            GD.debug("DIST: %s" % dist)
            return
        if dist <= 0.0:
            dist = 1.0
        self.camera.setDist(dist)
        self.camera.setClip(0.01*dist,100.*dist)
        self.camera.resetArea()

def structuredHexGrid(dx, dy, dz, isophex='hex64'):
    """it builds a structured hexahedral grid with nodes and elements both numbered in a structured way: first along z, then along y,and then along x. The resulting hex cells are oriented along z. This function is the equivalent of simple.rectangularGrid but for a mesh. Additionally, dx,dy,dz can be either integers or div (1D list or array). In case of list/array, first and last numbers should be 0.0 and 1.0 if the desired grid has to be inside the region 0.,0.,0. to 1.,1.,1.
    If isopHex is specified, a convenient set of control points for the isoparametric transformation hex64 is also returned.
    TODO: include other options to get the control points for other isoparametric transformation for hex."""
    sgx, sgy, sgz=dx, dy, dz
    if type(dx)!=int:sgx=len(dx)-1
    if type(dy)!=int:sgy=len(dy)-1
    if type(dz)!=int:sgz=len(dz)-1
    n3=regularGrid([0., 0., 0.],[1., 1., 1.],[sgx, sgy, sgz])
    if type(dx)!=int:n3[..., 0]=array(dx).reshape(-1, 1, 1)
    if type(dy)!=int:n3[..., 1]=array(dy).reshape(-1,  1)
    if type(dz)!=int:n3[..., 2]=array(dz).reshape(-1)
    nyz=(sgy+1)*(sgz+1)
    xh0= array([0, nyz, nyz+sgz+1,0+sgz+1 ])
    xh0= concatenate([xh0, xh0+1], axis=1)#first cell
    hz= array([xh0+j for j in range(sgz)])#z column
    hzy= array([hz+(sgz+1)*j for j in range(sgy)])#zy 2D rectangle
    hzyx=array([hzy+nyz*k for k in range(sgx)]).reshape(-1, 8)#zyx 3D
    if isophex=='hex64': return Coords(n3.reshape(-1, 3)), hzyx.reshape(-1, 8), regularGrid([0., 0., 0.], [1., 1., 1.], [3, 3, 3]).reshape(-1, 3)#control points for the hex64 applied to a basic struct hex grid
    else: return Coords(n3.reshape(-1, 3)), hzyx.reshape(-1, 8)

def run():
    reset()
    smooth()
    lights(True)

    S = TriSurface.read(getcfg('datadir')+'/horse.off')
    SA = draw(S)

    bb = S.bbox()
    bb1 = [ 1.1*bb[0]-0.1*bb[1], 1.1*bb[1]-0.1*bb[0]]
    print(bb)
    print(bb1)

    res = askItems([
        _I('Resolution',100),
        ])
    if not res:
        return
    
    nmax = res['Resolution']
    sz = bb1[1]-bb1[0]
    step = sz.max() / (nmax-1)
    n = (sz / step).astype(Int)
    print(n)
    P = Formex(simple.regularGrid(bb1[0],bb1[0]+n*step,n).reshape(-1,3))
    draw(P, marksize=1, color='black')
    #drawNumbers(P)
    zoomAll()
    ind = S.inside(P)
    vox = zeros(n+1,dtype=uint8)
    print(vox.shape)
    vox1 = vox.reshape(-1)
    print(vox1.shape,ind.max())
    vox1[ind] = 1
    print(vox.max())
    P.setProp(vox1)
    draw(P, marksize=8)

    dirname = askDirname()
    chdir(dirname)
    # Create output file
    if not checkWorkdir():
        print("Could not open a directory for writing. I have to stop here")
        return
    
    fs = utils.NameSequence('horse','.png')
    clear()
    flat()
    A = None
    for frame in vox:
        B = showGreyImage(frame)
        saveBinaryImage(frame*255,fs.next())
        undraw(A)
        A = B

def cone1(r0,r1,h,t=360.,nr=1,nt=24,diag=None):
    """Constructs a Formex which is (a sector of) a
    circle / (truncated) cone / cylinder.

    r0,r1,h are the lower and upper radius and the height of the truncated
    cone. All can be positive, negative or zero.
    Special cases:
    r0 = r1 : cylinder
    h = 0 : (flat) circle
    r0 = 0 or r1 = 0 : untruncated cone

    Only a sector of the structure, with opening angle t, is modeled.
    The default results in a full circumference.

    The cone is modeled by nr elements in height direction and nt elements in
    circumferential direction.
    
    By default, the result is a 4-plex Formex whose elements are quadrilaterals
    (some of which may collapse into triangles).
    If diag='up' or diag = 'down', all quads are divided by an up directed
    diagonal and a plex-3 Formex results.
    """
    r0,r1,h,t = map(float,(r0,r1,h,t))
    p = Formex(simple.regularGrid([r0,0.,0.],[r1,h,0.],[0,nr,0]).reshape(-1,3))
    #draw(p,color=red)
    a = (r1-r0)/h 
    if a != 0.:
        p = p.shear(0,1,a)
    #draw(p)
    q = p.rotate(t/nt,axis=1)
    #draw(q,color=green)
    if diag == 'u':
        F = connect([p,p,q],bias=[0,1,1]) + \
            connect([p,q,q],bias=[1,2,1])
    elif diag == 'd':
        F = connect([q,p,q],bias=[0,1,1]) + \
            connect([p,p,q],bias=[1,2,1])
    else:
        F = connect([p,p,q,q],bias=[0,1,1,0])

    F = Formex.concatenate([F.rotate(i*t/nt,1) for i in range(nt)])
    return F

    def setCamera(self,bbox=None,angles=None):
        """Sets the camera looking under angles at bbox.

        This function sets the camera parameters to view the specified
        bbox volume from the specified viewing direction.

        Parameters:

        - `bbox`: the bbox of the volume looked at
        - `angles`: the camera angles specifying the viewing direction.
          It can also be a string, the key of one of the predefined
          camera directions

        If no angles are specified, the viewing direction remains constant.
        The scene center (camera focus point), camera distance, fovy and
        clipping planes are adjusted to make the whole bbox viewed from the
        specified direction fit into the screen.

        If no bbox is specified, the following remain constant:
        the center of the scene, the camera distance, the lens opening
        and aspect ratio, the clipping planes. In other words the camera
        is moving on a spherical surface and keeps focusing on the same
        point.

        If both are specified, then first the scene center is set,
        then the camera angles, and finally the camera distance.

        In the current implementation, the lens fovy and aspect are not
        changed by this function. Zoom adjusting is performed solely by
        changing the camera distance.
        """
        #
        # TODO: we should add the rectangle (digital) zooming to
        #       the docstring

        self.makeCurrent()

        # set scene center
        if bbox is not None:
            pf.debug("SETTING BBOX: %s" % self.bbox,pf.DEBUG.DRAW)
            self.setBbox(bbox)

            X0,X1 = self.bbox
            self.camera.focus = 0.5*(X0+X1)

        # set camera angles
        if type(angles) is str:
            angles = self.view_angles.get(angles)
        if angles is not None:
            try:
                self.camera.setAngles(angles)
            except:
                raise ValueError,'Invalid view angles specified'

        # set camera distance and clipping planes
        if bbox is not None:
            # Currently, we keep the default fovy/aspect
            # and change the camera distance to focus
            fovy = self.camera.fovy
            #pf.debug("FOVY: %s" % fovy,pf.DEBUG.DRAW)
            self.camera.setLens(fovy,self.aspect)
            # Default correction is sqrt(3)
            correction = float(pf.cfg.get('gui/autozoomfactor',1.732))
            tf = coords.tand(fovy/2.)

            import simple
            bbix = simple.regularGrid(X0,X1,[1,1,1])
            bbix = dot(bbix,self.camera.rot[:3,:3])
            bbox = coords.Coords(bbix).bbox()
            dx,dy,dz = bbox[1] - bbox[0]
            vsize = max(dx/self.aspect,dy)
            dsize = bbox.dsize()
            offset = dz
            dist = (vsize/tf + offset) / correction

            if dist == nan or dist == inf:
                pf.debug("DIST: %s" % dist,pf.DEBUG.DRAW)
                return
            if dist <= 0.0:
                dist = 1.0
            self.camera.dist = dist

            ## print "vsize,dist = %s, %s" % (vsize,dist)
            ## near,far = 0.01*dist,100.*dist
            ## print "near,far = %s, %s" % (near,far)
            #near,far = dist-1.2*offset/correction,dist+1.2*offset/correction
            near,far = dist-1.0*dsize,dist+1.0*dsize
            # print "near,far = %s, %s" % (near,far)
            #print (0.0001*vsize,0.01*dist,near)
            # make sure near is positive
            near = max(near,0.0001*vsize,0.01*dist,finfo(coords.Float).tiny)
            # make sure far > near
            if far <= near:
                far += finfo(coords.Float).eps
            #print "near,far = %s, %s" % (near,far)
            self.camera.setClip(near,far)
            self.camera.resetArea()

"""

from plugins import isopar
import simple
import elements

wireframe()

ttype = ask("Select type of transformation",['Cancel','1D','2D','3D'])
if not ttype or ttype ==  'Cancel':
    exit()

tdim = int(ttype[0])

# create a unit quadratic grid in tdim dimensions
x = Coords(simple.regularGrid([0.]*tdim, [1.]*tdim, [2]*tdim)).reshape(-1,3)
x1 = Formex(x)
x2 = x1.copy()

# move a few points
if tdim == 1:
    eltype = 'line3'
    x2[1] = x2[1].rot(-22.5)
    x2[2] = x2[2].rot(22.5)
elif tdim == 2:
    eltype = 'quad9'
    x2[5] = x2[2].rot(-22.5)
    x2[8] = x2[2].rot(-45.)
    x2[7] = x2[2].rot(-67.5)
    x2[4] = x2[8] * 0.6
else:

                       (0,3,2,1,11,10,9,8), (4,5,6,7,12,13,14,15) ], ),
    reversed = (4,5,6,7,0,1,2,3,12,13,14,15,8,9,10,11,16,17,18,19),
)

Hex20.drawfaces = [ Hex20.faces.selectNodes(i) for i in Quad8.drawfaces ]
Hex20.drawfaces2 = [ Hex20.faces ]


# THIS ELEMENT USES A REGULAR NODE NUMBERING!!
# WE MIGHT SWITCH OTHER ELEMENTS TO THIS REGULAR SCHEME TOO
# AND ADD THE RENUMBERING TO THE FE OUTPUT MODULES
from simple import regularGrid
Hex27 = createElementType(
    'hex27',"A 27-node hexahedron",
    ndim = 3,
    vertices = regularGrid([0.,0.,0.],[1.,1.,1.],[2,2,2]).swapaxes(0,2).reshape(-1,3),
    edges = ('line3',[ (0,1,2),(6,7,8),(18,19,20),(24,25,26),
                       (0,3,6),(2,5,8),(18,21,24),(20,23,26),
                       (0,9,18),(2,11,20),(6,15,24),(8,17,26) ],),
    faces = ('quad9',[ (0,18,24,6,9,21,15,3,12),(2,8,26,20,5,17,23,11,14),
                       (0,2,20,18,1,11,19,9,10),(6,24,26,8,15,25,17,7,16),
                       (0,6,8,2,3,7,5,1,4),(18,20,26,24,19,23,25,21,22), ],),
)
Hex27.drawfaces = [ Hex27.faces.selectNodes(i) for i in Quad9.drawfaces ]

######################################################################
########## element type conversions ##################################

_dev_doc_ = """Element type conversion

Element type conversion in pyFormex is a powerful feature to transform