# I need to model an incline rectangle shape

**URL:** <https://fortran-lang.discourse.group/t/i-need-to-model-an-incline-rectangle-shape/3984>\
**Category:** Uncategorized\
**Created:** [July 12, 2022, 6:00am UTC](https://fortran-lang.discourse.group/t/i-need-to-model-an-incline-rectangle-shape/3984 "2022-07-12T06:00:25Z")\
**Posts on this page:** 1\
**Showing post:** 8

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [July 13, 2022, 1:12am UTC](https://fortran-lang.discourse.group/t/i-need-to-model-an-incline-rectangle-shape/3984/8 "2022-07-13T01:12:04Z")

</div>

Combining my two previous answers, omitting a few details, you can do something like this:

```fortran
   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:

```fortran
   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](http://www.computational-seismology.org/), 2016, by Heiner Igel.

 ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/e/e246a2d634e373a46060649fb1c692fb00104812.png)

---

_[View the full topic](https://fortran-lang.discourse.group/t/i-need-to-model-an-incline-rectangle-shape/3984)._