#! /usr/bin/env python # def tetrahedron_grid_points ( n, xv, ng ): #*****************************************************************************80 # ## TETRAHEDRON_GRID_POINTs computes points on a tetrahedral grid. # # Discussion: # # The grid is defined by specifying the coordinates of an enclosing # tetrahedron T, and the number of subintervals each edge of the # tetrahedron should be divided into. # # Choosing N = 10, for instance, breaks each side into 10 subintervals, # and produces a grid of ((10+1)*(10+2)*(10+3))/6 = 286 points. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 15 April 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the number of subintervals. # # Input, real XV[4,3], the vertices of the tetrahedron. # # Input, integer NG, the number of grid points. # # Output, real XG[NG,3], the tetrahedron grid points. # import numpy as np xg = np.zeros ( [ ng, 3 ] ) p = 0 for i in range ( 0, n + 1 ): for j in range ( 0, n + 1 - i ): for k in range ( 0, n + 1 - i - j ): l = n - i - j - k xg[p,0] = ( float ( i ) * xv[0,0] + float ( j ) * xv[1,0] \ + float ( k ) * xv[2,0] + float ( l ) * xv[3,0] ) / float ( n ) xg[p,1] = ( float ( i ) * xv[0,1] + float ( j ) * xv[1,1] \ + float ( k ) * xv[2,1] + float ( l ) * xv[3,1] ) / float ( n ) xg[p,2] = ( float ( i ) * xv[0,2] + float ( j ) * xv[1,2] \ + float ( k ) * xv[2,2] + float ( l ) * xv[3,2] ) / float ( n ) p = p + 1 return xg def tetrahedron_grid_points_test ( ): #*****************************************************************************80 # ## TETRAHEDRON_GRID_POINTS_TEST tests TETRAHEDRON_GRID_POINTS. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 10 November 2011 # # Author: # # John Burkardt # import numpy as np import platform from r83col_print_part import r83col_print_part from r8mat_write import r8mat_write from tetrahedron_grid_count import tetrahedron_grid_count from tetrahedron_grid_display import tetrahedron_grid_display print ( '' ) print ( 'TETRAHEDRON_GRID_POINTS_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' TETRAHEDRON_GRID_POINTS can define a tetrahedral grid' ) print ( ' with N+1 points on a side, based on any tetrahedron.' ) n = 5 print ( ' N = %d' % ( n ) ) ng = tetrahedron_grid_count ( n ) # # Define the tetrahedron vertices. # xv = np.array ( [ \ [ 0.0, 0.0, 0.0 ], \ [ 2.0, 1.0, 0.0 ], \ [ 1.0, 4.0, 0.0 ], \ [ 3.0, 3.0, 3.0 ] ] ) r83col_print_part ( 4, xv, 20, ' Tetrahedron vertices:' ) xg = tetrahedron_grid_points ( n, xv, ng ) r83col_print_part ( ng, xg, 20, ' Tetrahedron grid points:' ) # # Write the data to a file. # filename = 'tetrahedron_grid_points.xyz' r8mat_write ( filename, ng, 3, xg ) print ( '' ) print ( ' Data written to the file "%s".' % ( filename ) ) # # Plot the data. # filename = 'tetrahedron_grid_points.png' tetrahedron_grid_display ( xv, ng, xg, filename ) print ( '' ) print ( ' Plot written to the file "%s".' % ( filename ) ) # # Terminate. # print ( '' ) print ( 'TETRAHEDRON_GRID_POINTS_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) tetrahedron_grid_points_test ( ) timestamp ( )