diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index 552b00bd5fe4aafa737fedd9364023585b933c64..ae35a62455b680e19f874613b2b2cda63e5bf1d3 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -16,7 +16,7 @@ PUBLIC :: init_ellip, boundary_conditions_ellip, postprocess_ellip, &
contains
-subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
+subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,dz,remp)
use decomp_2d, only : mytype
use param, only : one, two, ten
@@ -30,14 +30,14 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
integer :: nxi,nxf,ny,nyi,nyf,nzi,nzf
real(mytype),dimension(nxi:nxf,nyi:nyf,nzi:nzf) :: epsi
real(mytype),dimension(ny) :: yp
- real(mytype) :: dx
+ real(mytype) :: dx,dz
real(mytype) :: remp
integer :: i,j,k
real(mytype) :: xm,ym,zm,r,rads2,kcon
real(mytype) :: zeromach
real(mytype) :: cexx,ceyy,cezz,dist_axi,eqr
real(mytype) :: point(3)
- logical :: is_inside, is_inside_2
+ logical :: is_inside, is_inside_2, is_inside_3
zeromach=one
do while ((one + zeromach / two) .gt. one)
@@ -67,6 +67,18 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
linearVelocity_2=[lvx2,lvy2,lvz2]
angularVelocity_2=[zero,avx2,avy2,avz2]
call ellipInertiaCalculate(shape_2,rho_s_2,inertia_2)
+ if (nobjmax.gt.2) then
+ eqr=(shx3*shy3*shz3)**(1.0/3.0)
+ shape_3=[shx3/eqr,shy3/eqr,shz3/eqr]
+
+ orientation_3=[oriw3,orii3,orij3,orik3]
+ call NormalizeQuaternion(orientation_3)
+ ! write(*,*) 'Quaternion normalized!'
+ position_3=[cex3,cey3,cez3]
+ linearVelocity_3=[lvx3,lvy3,lvz3]
+ angularVelocity_3=[zero,avx3,avy3,avz3]
+ call ellipInertiaCalculate(shape_3,rho_s_3,inertia_3)
+ endif
endif
endif
@@ -110,7 +122,10 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
call is_inside_ellipsoid(point, position, orientation, shape, ra, zeromach, is_inside)
if (nobjmax.gt.1) then
- call is_inside_ellipsoid(point, position_2, orientation_2, shape_2, ra, zeromach, is_inside_2)
+ call is_inside_ellipsoid(point, position_2, orientation_2, shape_2, ra2, zeromach, is_inside_2)
+ if (nobjmax.gt.2) then
+ call is_inside_ellipsoid(point, position_3, orientation_3, shape_3, ra3, zeromach, is_inside_3)
+ endif
endif
! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two+(zm-cezz)**two)
@@ -123,7 +138,9 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
if ((.not.is_inside).and.(nobjmax.eq.1)) then
! write(*,*) i, j, k
cycle
- else if ((.not.is_inside).and.(.not.is_inside_2)) then
+ else if ((.not.is_inside).and.(.not.is_inside_2).and.(nobjmax.eq.2)) then
+ cycle
+ else if ((.not.is_inside).and.(.not.is_inside_2).and.(.not.is_inside_3).and.(nobjmax.eq.3)) then
cycle
endif
! write(*,*) i, j, k, zm
diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90
index f53b8a6a9bd5d57c7e92edb37a7e0b023f88bcb6..416231f95b8024f3d70589330ae9976f54cfb11f 100644
--- a/src/genepsi3d.f90
+++ b/src/genepsi3d.f90
@@ -79,7 +79,7 @@ contains
ELSEIF (itype.EQ.itype_ellip) THEN
- CALL geomcomplex_ellip(epsi, nxi, nxf, ny, nyi, nyf, nzi, nzf, dx, yp, remp)
+ CALL geomcomplex_ellip(epsi, nxi, nxf, ny, nyi, nyf, nzi, nzf, dx, yp, dz, remp)
ENDIF
diff --git a/src/module_param.f90 b/src/module_param.f90
index 63466ccb75d4e38807369565e7179d4acd85013b..a2c814b2362f41255cea8f3b03ad8182ccabf694 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -618,11 +618,15 @@ end module simulation_stats
module ibm_param
use decomp_2d, only : mytype
real(mytype) :: cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,ubcx,ubcy,ubcz,rads,rho_s,c_air
- real(mytype) :: cex2,cey2,cez2,shx2,shy2,shz2,oriw2,orii2,orij2,orik2,lvx2,lvy2,lvz2,avx2,avy2,avz2,rho_s_2
+ real(mytype) :: cex2,cey2,cez2,shx2,shy2,shz2,oriw2,orii2,orij2,orik2,lvx2,lvy2,lvz2,avx2,avy2,avz2,rho_s_2,ra2
+ real(mytype) :: cex3,cey3,cez3,shx3,shy3,shz3,oriw3,orii3,orij3,orik3,lvx3,lvy3,lvz3,avx3,avy3,avz3,rho_s_3,ra3
real(mytype) :: position(3),orientation(4),linearVelocity(3),angularVelocity(4),linearAcceleration(3),torque(3),shape(3),inertia(3,3)
real(mytype) :: position_2(3),orientation_2(4),linearVelocity_2(3),angularVelocity_2(4),linearAcceleration_2(3),torque_2(3),shape_2(3),inertia_2(3,3)
+ real(mytype) :: position_3(3),orientation_3(4),linearVelocity_3(3),angularVelocity_3(4),linearAcceleration_3(3),torque_3(3),shape_3(3),inertia_3(3,3)
+
real(mytype) :: position_1(3),linearVelocity_1(3),orientation_1(4),angularVelocity_1(4)
real(mytype) :: position_1_2(3),linearVelocity_1_2(3),orientation_1_2(4),angularVelocity_1_2(4)
+ real(mytype) :: position_1_3(3),linearVelocity_1_3(3),orientation_1_3(4),angularVelocity_1_3(4)
real(mytype) :: chord,thickness,omega
integer :: inana ! Analytical BC as Input
integer :: imove
diff --git a/src/navier.f90 b/src/navier.f90
index d18cc00655cf6eee300b3df50dbd621cece6c3c4..e01003dee97c12d1b8075ec36bfc461b1996ac1f 100644
--- a/src/navier.f90
+++ b/src/navier.f90
@@ -280,6 +280,8 @@ contains
real(mytype),dimension(ph1%zst(1):ph1%zen(1),ph1%zst(2):ph1%zen(2),nzmsize) :: pp3
real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1_ux,ep1_uy,ep1_uz
real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1_ux2,ep1_uy2,ep1_uz2
+ real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1_ux3,ep1_uy3,ep1_uz3
+
integer :: nvect3,i,j,k,nlock
integer :: code
@@ -304,6 +306,15 @@ contains
ta1(:,:,:) = (one - ep1(:,:,:)) * ux1(:,:,:) + ep1(:,:,:)*(ep1_ux(:,:,:)+ep1_ux2(:,:,:))
tb1(:,:,:) = (one - ep1(:,:,:)) * uy1(:,:,:) + ep1(:,:,:)*(ep1_uy(:,:,:)+ep1_uy2(:,:,:))
tc1(:,:,:) = (one - ep1(:,:,:)) * uz1(:,:,:) + ep1(:,:,:)*(ep1_uz(:,:,:)+ep1_uz2(:,:,:))
+
+ else if (nobjmax.eq.3) then
+ call navierFieldGen(position, linearVelocity, angularVelocity,ep1, ep1_ux, ep1_uy, ep1_uz)
+ call navierFieldGen(position_2, linearVelocity_2, angularVelocity_2, ep1, ep1_ux2, ep1_uy2, ep1_uz2)
+ call navierFieldGen(position_3, linearVelocity_3, angularVelocity_3, ep1, ep1_ux3, ep1_uy3, ep1_uz3)
+
+ ta1(:,:,:) = (one - ep1(:,:,:)) * ux1(:,:,:) + ep1(:,:,:)*(ep1_ux(:,:,:)+ep1_ux2(:,:,:)+ep1_ux3(:,:,:))
+ tb1(:,:,:) = (one - ep1(:,:,:)) * uy1(:,:,:) + ep1(:,:,:)*(ep1_uy(:,:,:)+ep1_uy2(:,:,:)+ep1_uy3(:,:,:))
+ tc1(:,:,:) = (one - ep1(:,:,:)) * uz1(:,:,:) + ep1(:,:,:)*(ep1_uz(:,:,:)+ep1_uz2(:,:,:)+ep1_uz3(:,:,:))
else
write(*,*) 'Currently unsupported number of bodies'
endif
diff --git a/src/parameters.f90 b/src/parameters.f90
index 9bcddb98ecd96ed9c20a3f5f769076d35ab93c24..e0dddf1d85cb1591c84066dd8a3b2771832ed632 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -61,7 +61,8 @@ subroutine parameter(input_i3d)
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,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,nobjmax, &
- cex2,cey2,cez2,shx2,shy2,shz2,oriw2,orii2,orij2,orik2,lvx2,lvy2,lvz2,avx2,avy2,avz2,rho_s_2, &
+ cex2,cey2,cez2,shx2,shy2,shz2,oriw2,orii2,orij2,orik2,lvx2,lvy2,lvz2,avx2,avy2,avz2,rho_s_2,ra2, &
+ cex3,cey3,cez3,shx3,shy3,shz3,oriw3,orii3,orij3,orik3,lvx3,lvy3,lvz3,avx3,avy3,avz3,rho_s_3,ra3, &
nraf,nvol,iforces, npif, izap, ianal, imove, thickness, chord, omega ,ubcx,ubcy,ubcz,rads,rho_s, c_air
NAMELIST /ForceCVs/ xld, xrd, yld, yud!, zld, zrd
NAMELIST /LMN/ dens1, dens2, prandtl, ilmn_bound, ivarcoeff, ilmn_solve_temp, &
@@ -573,7 +574,7 @@ subroutine parameter_defaults()
use probes, only : nprobes, flag_all_digits, flag_extra_probes
use visu, only : output2D
use forces, only : iforces, nvol
- use ibm_param, only : oriw,orii,orij,orik,rho_s,shx,shy,shz,oriw2,shx2,shy2,shz2,rho_s_2
+ use ibm_param, only : oriw,orii,orij,orik,rho_s,shx,shy,shz,oriw2,shx2,shy2,shz2,rho_s_2,oriw3,shx3,shy3,shz3,rho_s_3
implicit none
@@ -624,8 +625,13 @@ subroutine parameter_defaults()
shx2=one
shy2=one
shz2=one
+ shx3=one
+ shy3=one
+ shz3=one
oriw2=one
+ oriw3=one
rho_s_2=one
+ rho_s_3=one
!! Gravity field
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index 8f5e98eeccb70d28f7efd6845b14991002da5e3f..bab91245c1e2cda131e21ef0d2e56c76d9a85f8f 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -104,6 +104,20 @@ program xcompact3d
orientation_2 = orientation_1_2
angularVelocity_2 = angularVelocity_1_2
+ if (nobjmax.gt.2) then
+ linearAcceleration_3(:)=zero
+ torque_3(:)=zero
+
+ call lin_step(position_3,linearVelocity_3,linearAcceleration_3,dt,position_1_3,linearVelocity_1_3)
+ call ang_step(orientation_3,angularVelocity_3,torque_3,dt,orientation_1_3,angularVelocity_1_3)
+
+ position_3 = position_1_3
+ linearVelocity_3 = linearVelocity_1_3
+
+ orientation_3 = orientation_1_3
+ angularVelocity_3 = angularVelocity_1_3
+
+ end if
end if
endif