nx=1001
nz=501
idepth=30
allocate(icode(1:nx,1:nz))
do i=1,nx
do j=1,nz
if(j<=idepth)then
icode(i,j)='1'
elseif(j>idepth)then
icode(i,j)='2'
endif
enddo
enddo
do i=1,nx
do j=1,nz
if(i<=525.and.i>=475.and.j<=14.and.j>=5)then
icode(i,j)='3'
endif
enddo
enddo
This is starting to look like a homework problem. The parallelogram has four sides. The top and bottom are simple conditions. The sides are each determined by a linear equation. A point is inside if it satisfies four test conditions. The order of the tests does not matter for the result, but it does matter as far as efficiency. You want the failures to be recognized as early as possible, so arrange the four tests in that way. This depends on how the population of points to be tested is distributed in the overall grid.
Have you considered investing some time to study computational geometry? Your program would be much more flexible than now when you are manually encoding tests for each geometry variation. It could save you a lot of compilation and debugging.
Say you take the PNPOLY routine by W. Randolph Franklin. Compile it separately from the rest of your code with the command gfortran -std=legacy -c pnpoly.f.
Then you could perform your test as follows:
interface
subroutine pnpoly(px,py,xx,yy,n,inout)
implicit none
integer, intent(in) :: n
real, intent(in) :: px, py
real, intent(in) :: xx(n), yy(n)
integer, intent(out) :: inout
end subroutine
end interface
integer, parameter :: wp = kind(1.0e0)
integer, parameter :: nx = 101, nz = 51
real(wp) :: hx, hz, xx, zz
integer :: poly_n, in_polygon
real(wp), allocatable :: poly_x(:), poly_z(:)
! step size
hz = 0.5_wp
hx = 1.0_wp
! describe polygon as x- and z-coordinates of vertices
! (you could read these from a namelist or some geometry file)
poly_n = 4
poly_x = [ ... ]
poly_z = [ ... ]
! for each grid point...
do i = 1, nx
do j = 1, nz
! transform from grid to spatial coordinates
xx = (i-1)*hx
zz = 0 - (j-1)*hz
! perform point-in polygon test
call pnpoly(xx,zz,poly_x,poly_z,poly_n,in_polygon)
if (in_polygon >= 0) then
! inside polygon
icode(i,j) = '3'
end if
end do
end do
The question said âinclined rectangleâ but the picture showed a parallelogram. To me an inclined rectangle would have its sides at some angle not a multiple of 90 degrees to the vertical and horizontal but at 0 or 90 degrees to one another.
Well the point in polygon approach should work regardless of the type of quadrilateral.
I pulled a few bits and pieces from a few programs to create VTK output which can be visualized in Paraview. Consider this only as an untested prototype which should be adapted further to your needs. I was using the XY-coordinates, but to maintain notation cleanliness probably X and Z should be used. Also I wasnât sure what reference frame you have for spatial coordinates.
module precision
implicit none
public
!> Working-precision real kind
integer, parameter :: wp = kind(1.0e0)
end module
module vtk_format
use precision, only: wp
implicit none
public :: FMT_REAL, FMT_VECTOR
!> Format string for single precision real numbers
character(*), parameter :: FMT_SP = '(ES15.8E2)'
!> Format string for double precision real numbers
character(*), parameter :: FMT_DP = '(ES24.16E3)'
!> Format string for working-precision real kind
character(*), parameter :: FMT_REAL = &
trim(merge(FMT_SP//' ', FMT_DP, wp == kind(1.0e0)))
!> Format string for a vector triplet
character(*), parameter :: FMT_VECTOR = '(3'//FMT_REAL//')'
end module
module grid
use precision, only: wp
use vtk_format, only: FMT_REAL, FMT_VECTOR
implicit none
private
public :: output_vtk_structuredPoints
contains
subroutine output_vtk_structuredPoints(filename,nx,ny,hx,hy,rho,description)
character(len=*), intent(in) :: filename
!> Output file name
integer, intent(in) :: nx, ny
!> Mesh dimensions, number of cells in each direction
!> the number of grid points is nx + 1, ny + 1
real(wp), intent(in) :: hx, hy
!> Grid spacing, must be larger than 0
real(wp), intent(in) :: rho(:,:)
!> A scalar field defined at the grid points
character(len=*), intent(in), optional :: description
!> Optionally, single-line description of the data
integer :: unit, i, j
open(newunit=unit,file=filename)
write(unit, '(A)') "# vtk DataFile Version 3.0"
if (present(description)) then
write(unit, '(A)') description
else
write(unit, '(A)') "scalar_field"
end if
write(unit, '(A)') "ASCII"
write(unit, '(A)') "DATASET STRUCTURED_POINTS"
write(unit, '(A,3(I0,:,1X))') "DIMENSIONS ", nx + 1, ny + 1, 1
write(unit, '(A,'// FMT_VECTOR //')') "ORIGIN ", 0.0_wp, 0.0_wp, 0.0_wp
write(unit, '(A,'// FMT_VECTOR //')') "SPACING ", &
hx, hy, min(hx, hy)
!
! Dataset attributes
!
write(unit, *)
write(unit, '(A,I0)') "POINT_DATA ", (nx + 1)*(ny + 1)*1
write(unit, '(A)') "SCALARS rho float 1"
write(unit, '(A)') "LOOKUP_TABLE default"
do j = 1, ny + 1
do i = 1, nx + 1
write(unit, FMT_REAL) rho(i,j)
end do
end do
close(unit)
end subroutine output_vtk_structuredPoints
end module
module polygons
use precision, only: wp
use vtk_format, only: FMT_REAL, FMT_VECTOR
implicit none
private
public :: point, polygon, polygon_to_vtk
type :: point
real(wp) :: x, y
end type
type :: polygon
integer :: nvert
real(wp), allocatable :: x(:), y(:)
end type
interface polygon
module procedure :: new_polygon
end interface
contains
function new_polygon(points) result(res)
type(point) :: points(:)
type(polygon) :: res
integer :: i
res%nvert = size(points)
associate(n => res%nvert)
allocate(res%x(n),res%y(n))
do i = 1, n
res%x(i) = points(i)%x
res%y(i) = points(i)%y
end do
end associate
end function
subroutine polygon_to_vtk(filename,poly)
character(len=*), intent(in) :: filename
type(polygon), intent(in) :: poly
integer, parameter :: VTK_POLYGON = 7
integer :: unit, i
open(newunit=unit,file=filename)
write(unit, '(A)') "# vtk DataFile Version 3.0"
write(unit, '(A)') ""
write(unit, '(A)') "ASCII"
write(unit, '(A)') "DATASET UNSTRUCTURED_GRID"
write(unit, '(A, I0, A)') "POINTS ", poly%nvert, " float"
do i = 1, poly%nvert
write(unit, '('// FMT_VECTOR //')') poly%x(i), poly%y(i), 0.0_wp
end do
write(unit, '(A, 2(I0,:,2X))') "CELLS ", 1, 1+poly%nvert
write(unit, '(*(I0,:,1X))') poly%nvert, (i-1, i = 1, poly%nvert)
write(unit, '(A, I0)') "CELL_TYPES ", 1
write(unit, '(I0)') VTK_POLYGON
end subroutine
end module
program test
use precision, only: wp
use grid, only: output_vtk_structuredPoints
use polygons
implicit none
integer, parameter :: nx = 101, nz = 51
real(wp) :: hx, hz
real(wp), allocatable :: rho(:,:)
type(polygon) :: obstacle
hx = 1.0_wp
hz = 1.0_wp
allocate(rho(nx,nz))
rho = 1.0_wp
call output_vtk_structuredPoints("test.vtk",nx-1,nz-1,hx,hz,rho,"rho")
obstacle = polygon([point(40._wp, 10._wp), &
point(60._wp, 10._wp), &
point(70._wp, 30._wp), &
point(50._wp, 30._wp)])
call to_vtk("obstacle.vtk",obstacle)
end program
Thank you very much
Combining my two previous answers, omitting a few details, you can do something like this:
impure function point_in_polygon(px,py,poly) result(res)
real(wp), intent(in) :: px, py
type(polygon), intent(in) :: poly
interface
subroutine pnpoly(px,py,xx,yy,n,inout)
implicit none
integer, intent(in) :: n
real, intent(in) :: px, py
real, intent(in) :: xx(n), yy(n)
integer, intent(out) :: inout
end subroutine
end interface
logical :: res
integer :: ires
res = .false.
call pnpoly(px,py,poly%x,poly%y,poly%nvert,ires)
if (ires >= 0) then
res = .true.
end if
end function
In the main program (or preferably a subroutine) youâd apply the filter:
rho = 1.0_wp
do j = 1, nz
do i = 1, nx
xx = (i-1)*hx
zz = (j-1)*hz
if (point_in_polygon(xx,zz,obstacle)) then
rho(i,j) = 2.0_wp
end if
end do
end do
call output_vtk_structuredPoints("test.vtk",nx-1,nz-1,hx,hz,rho,"rho")
After refreshing Paraview you can see the point-in-polygon works, however there is some ambiguity at the edges. A pragmatic solution would be to apply a âfudgeâ factor to the sides of the quadrilateral. This should be a value like 10*epsilon(xx). It might also be an error in the PNPOLY routine. But generally speaking itâs a problem of your discretization method. If you are simulating wave propagation in bodies with different properties, you are bound to commit some errors at the interface. You can find a few other approaches in the excellent book, Computational Seismology: A Practical Introduction, 2016, by Heiner Igel.
I want to create an oblique rectangle
implicit none
integer:: i,j,nx,nz
integer :: idepth
character*1,allocatable,dimension(:,
:: icode
nx=1001
nz=501
idepth=30
allocate(icode(1:nx,1:nz))
do i=1,nx
do j=1,nz
if(j<=idepth)then
icode(i,j)='1'
elseif(j>idepth)then
icode(i,j)='2'
endif
enddo
enddo
do i=1,nx
do j=1,nz
if(i<=525.and.i>=475.and.j<=14.and.j>=5)then
icode(i,j)='3'
endif
enddo
enddo
open(2012,file=âmodel.txtâ,status=âunknownâ)
do j=1,nz
write(2012,1000)(icode(i,j),i=1,nx)
enddo
close(2012)
1000 format(a1)
How to modify this one? This is a normal rectangular shape. I need to change it to be inclined
Use the general form of a line a*x + b*y + c = 0 (opposed to the more conventional slope-intercept form). The line divides the plane into two parts, the part where a*x + b*y + c > 0, and the part with a*x + b*y + c < 0.
Two of your lines are parallel to the x-axis, so those are easy (y = L3, and y = L4). For the other two, you need to find the coefficients (a,b,c).
implicit none
integer:: i,j,nx,nz
integer:: a,b,c,x,y
integer:: L1,L2,L3,L4
integer :: idepth
character*1,allocatable,dimension(:, :: icode
nx=1001
nz=501
idepth=30
x=500
y=100
L1=400
L2=520
L3=520
L4=600
a= 30
b= 20
c= 20
allocate(icode(1:nx,1:nz))
do i=1,nx
do j=1,nz
if(j<=idepth)then
icode(i,j)='1'
elseif(j>idepth)then
icode(i,j)='2'
endif
enddo
enddo
do i=1,nx
do j=1,nz
if(((a*x) + (b*y )+ (c = 0)).and. ((a*x) + (b*y) + (c > 0)).and. ((a*x) + (b*y) + (c < 0)))then
icode(i,j)='3'
endif
enddo
enddo
open(2012,file=âmodel.txtâ,status=âunknownâ)
do j=1,nz
write(2012,1000)(icode(i,j),i=1,nx)
enddo
close(2012)
1000 format(a1)