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.
