From f9f04631c96d583ad4e546e2945ccc7ea7718c8a Mon Sep 17 00:00:00 2001
From: Andrew Boyd <adfboyd@gmail.com>
Date: Tue, 14 Feb 2023 13:55:50 +0000
Subject: [PATCH] added variables for elliptical shape

---
 src/BC-Cylinder.f90  | 15 ++++++++++-----
 src/module_param.f90 |  2 +-
 src/parameters.f90   |  2 +-
 3 files changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/BC-Cylinder.f90 b/src/BC-Cylinder.f90
index 749a73f3..20eb4466 100644
--- a/src/BC-Cylinder.f90
+++ b/src/BC-Cylinder.f90
@@ -35,9 +35,9 @@ contains
     real(mytype)               :: dx
     real(mytype)               :: remp
     integer                    :: i,j,k
-    real(mytype)               :: xm,ym,r,rads2,kcon
+    real(mytype)               :: xm,ym,zm,r,rads2,kcon
     real(mytype)               :: zeromach
-    real(mytype)               :: cexx,ceyy,dist_axi
+    real(mytype)               :: cexx,ceyy,cezz,dist_axi
 
     zeromach=one
     do while ((one + zeromach / two) .gt. one)
@@ -48,30 +48,35 @@ contains
     ! Intitialise epsi
     epsi(:,:,:)=zero
 
-    ! Update center of moving Cylinder
+    ! Update center of moving Cylinder - change to ellipsoid
     !cexx=cex+ubcx*t
     !ceyy=cey+ubcy*t
     ! Update center of moving Cylinder
     if (t.ne.0.) then
        cexx=cex+ubcx*(t-ifirst*dt)
        ceyy=cey+ubcy*(t-ifirst*dt)
+       cezz=cez+ubcz*(t-ifirst*dt)
     else
        cexx=cex
        ceyy=cey
+       cezz=cez
     endif
     !
     ! Define adjusted smoothing constant
 !    kcon = log((one-0.0001)/0.0001)/(smoopar*0.5*dx) ! 0.0001 is the y-value, smoopar: desired number of affected points 
 !
     do k=nzi,nzf
+       zm=real(k+xstart(3)-2,mytype)*dz
        do j=nyi,nyf
           ym=yp(j)
           do i=nxi,nxf
              xm=real(i-1,mytype)*dx
-             r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two)
+             r=sqrt_prec(((xm-cexx)/sex)**two+((ym-ceyy)/sey)**two+((zm-cezz)/sez)**two)
              if (r-ra.gt.zeromach) then
-                cycle
+                
+               cycle
              endif
+             write(*,'(A,3(F5.2,2x))') "Inside solid, the coords are: ", xm, ym, zm
              epsi(i,j,k)=remp
           enddo
        enddo
diff --git a/src/module_param.f90 b/src/module_param.f90
index 49b0303c..85cb4ac9 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -616,7 +616,7 @@ end module simulation_stats
 !############################################################################
 module ibm_param
   use decomp_2d, only : mytype
-  real(mytype) :: cex,cey,cez,ra,ubcx,ubcy,ubcz,rads, c_air
+  real(mytype) :: cex,cey,cez,ra,sex,sey,sez,ubcx,ubcy,ubcz,rads, c_air
   real(mytype) :: chord,thickness,omega
   integer :: inana ! Analytical BC as Input
   integer :: imove
diff --git a/src/parameters.f90 b/src/parameters.f90
index 74b19756..0a551f65 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -60,7 +60,7 @@ subroutine parameter(input_i3d)
   NAMELIST /LESModel/ jles, smagcst, smagwalldamp, nSmag, walecst, maxdsmagcst, iwall
   NAMELIST /WallModel/ smagwalldamp
   NAMELIST /Tripping/ itrip,A_tr,xs_tr_tbl,ys_tr_tbl,ts_tr_tbl,x0_tr_tbl
-  NAMELIST /ibmstuff/ cex,cey,cez,ra,nobjmax,nraf,nvol,iforces, npif, izap, ianal, imove, thickness, chord, omega ,ubcx,ubcy,ubcz,rads, c_air
+  NAMELIST /ibmstuff/ cex,cey,cez,ra,sex,sey,sez,nobjmax,nraf,nvol,iforces, npif, izap, ianal, imove, thickness, chord, omega ,ubcx,ubcy,ubcz,rads, c_air
   NAMELIST /ForceCVs/ xld, xrd, yld, yud!, zld, zrd
   NAMELIST /LMN/ dens1, dens2, prandtl, ilmn_bound, ivarcoeff, ilmn_solve_temp, &
        massfrac, mol_weight, imultispecies, primary_species, &
-- 
GitLab