**Meshing information**
         Pts        SubGrid     Neighbor    Internal Edge    OutFlow Edge
         0          0(3,4,0)     (1,3)        0(3,4)            0(3,0)
         1          1(0,4,1)     (2,0)        1(4.0)            1(0,1)
         2          2(1,4.2)     (3,1)        2(1,4)            2(1,2)
         3          3(2,4,3)     (0,2)        3(2,4)            3(2,3)
         4

#1. Initialize a Wellblock with subdivision method
    subdivision.WellGrid(Rect0=(0,0),Rect1=(50,50),Qw=1000,Qe=(200,400,300,100),h=26.25,phi=0.2)
"""
        Rect0  -- The coordinates of left bottom corner point of the rectangle
        Rect1  -- The coordinates of right top corner point of the rectangle
        Qw     -- Wellbore flow rate, stb/day
        Qe     -- Outflow edge flow rate, stb/day (Edge0,Edge1,Edge2,Edge3)
        h      -- Wellblock thickness
        phi    -- Wellblock porosity
"""
#2. Solve the internal mass balnaced equation
    subdivision.Subdivision()
#3. Tracing streamline and output streamline and TOF information
    SL,TOF,SL_end,TOF_end=subdivision.SLTrace(NSL=80)
"""
        NSL        -- Number of Streamline
        TOF        -- TOF array of of the [NSL] streamlines  
        SL         -- nodes array of of the [NSL] streamlines  
        TOF_end    -- a boundary TOF list of the [NSL] streamlines  
        SL_end     -- a boundary node list of the [NSL] streamlines 
"""

**Meshing information**
         Pts        Outflow Edge    Well    
         0          0(0,1)        0
         1          1(1,2)        1
         2          2(2,3)        ..
         3          3(3,0)
         ..         ..

#1. Initialize a Wellblock with the embedded method
    embedded.WellGrid(Pts_e=Boundary_vert,Pts_w=Well_vert,Qe=Qedge,Qw=Qwell,Nbd=10,rw=0.25,h=26.25,phi=0.2,miu=1,kxy=(200,200))
"""
        Pts_e  -- Vertices of boundary edge e.g. [(0,0),(0,1),(1,1),(1,0)]
        Pts_w  -- Location of well e.g. [(0.5,0.5),(0.75,0.75)]
        Qw     -- Wellbore flow rate, stb/day e.g. [1000,200]
        Qe     -- Outflow edge flow rate, stb/day e.g. [400,300,300,200]
        h      -- Wellblock thickness
        phi    -- Wellblock porosity
        miu    -- Fluid viscosity
        kxy    -- Principle permeability of the wellblock. e.g. [kxx,kyy]    
"""
#2. Construct mesh for VBEM
    embedded.Meshing()
#3. Solve the velocity field for VBEM
    embedded.FlowSol()
#4. Tracing streamline and output streamline and TOF information
    SL,TOF,SL_end,TOF_end=embedded.SLtrace(NSL=80,deltaT=0.1,method='RK4')
"""
        NSL        -- Number of Streamline
        method      -- Numerical integration method. RK2, RK4 and Adapative are provided
        TOF        -- TOF array of of the [NSL] streamlines  
        SL         -- nodes array of of the [NSL] streamlines  
        TOF_end    -- a boundary TOF list of the [NSL] streamlines  
        SL_end     -- a boundary node list of the [NSL] streamlines 
"""
#6. Plotting and output the velocity field
    P,V_x,V_y=embedded.FieldPlot(vmax=50)
"""
        vmax     -- Maximum range of the color bar
        P        -- pressure array of the polygon domain  
        V_x      -- velocity in x direction of the polygon domain   
        V_y      -- velocity in y direction of the polygon domain  
"""