@@ -258,34 +258,34 @@ contains
258258
259259                    ! Calculate velocity of ghost cell
260260                    if  (gp%slip) then 
261-                       norm(1 :3 ) =  levelset_norm%sf(gp%loc(1 ), gp%loc(2 ), gp%loc(3 ), gp%ib_patch_id, 1 :3 )
262-                       buf =  sqrt (sum (norm** 2 ))
263-                       norm =  norm/ buf
264-                       vel_norm_IP =  sum (vel_IP* norm)* norm
265-                       vel_g =  vel_IP -  vel_norm_IP
266-                       if  (patch_ib(patch_id)%moving_ibm /=  0 ) then 
267-                           ! compute the linear velocity of the ghost point due to  rotation
268-                           radial_vector =  physical_loc -  [patch_ib(patch_id)%x_centroid, &
269-                                                           patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid]
270-                           rotation_velocity =  cross_product(matmul (patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector)
271- 
272-                           ! add only the component of the IB' s motion that is normal to the surface
273-                           vel_g = vel_g + sum((patch_ib(patch_id)%vel + rotation_velocity)*norm)*norm 
274-                       end if 
275-                   else 
276-                       if (patch_ib(patch_id)%moving_ibm == 0) then 
277-                           ! we know the object is not moving if moving_ibm is 0 (false) 
278-                           vel_g = 0._wp 
279-                       else 
280-                           ! get the vector that points from the centroid to the ghost 
281-                           radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, & 
282-                                                           patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid] 
283-                           ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear velocity 
284-                           rotation_velocity = cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector) 
285-                           do q = 1, 3 
286-                               ! if mibm is 1 or 2, then the boundary may be moving 
287-                               vel_g(q) = patch_ib(patch_id)%vel(q) ! add the linear velocity 
288-                               vel_g(q) = vel_g(q) + rotation_velocity(q) ! add the rotational velocity 
261+                          norm(1 :3 ) =  levelset_norm%sf(gp%loc(1 ), gp%loc(2 ), gp%loc(3 ), gp%ib_patch_id, 1 :3 )
262+                          buf =  sqrt (sum (norm** 2 ))
263+                          norm =  norm/ buf
264+                          vel_norm_IP =  sum (vel_IP* norm)* norm
265+                          vel_g =  vel_IP -  vel_norm_IP
266+                          if  (patch_ib(patch_id)%moving_ibm /=  0 ) then 
267+                              ! compute the linear velocity of the ghost point due to  rotation
268+                              radial_vector =  physical_loc -  [patch_ib(patch_id)%x_centroid, &
269+                                                              patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid]
270+                              rotation_velocity =  cross_product(matmul (patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector)
271+ 
272+                              ! add only the component of the IB' s motion that is normal to the surface
273+                              vel_g = vel_g + sum((patch_ib(patch_id)%vel + rotation_velocity)*norm)*norm 
274+                          end if 
275+                      else 
276+                          if (patch_ib(patch_id)%moving_ibm == 0) then 
277+                              ! we know the object is not moving if moving_ibm is 0 (false) 
278+                              vel_g = 0._wp 
279+                          else 
280+                              ! get the vector that points from the centroid to the ghost 
281+                              radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, & 
282+                                                              patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid] 
283+                              ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear velocity 
284+                              rotation_velocity = cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector) 
285+                              do q = 1, 3 
286+                                  ! if mibm is 1 or 2, then the boundary may be moving 
287+                                  vel_g(q) = patch_ib(patch_id)%vel(q) ! add the linear velocity 
288+                                  vel_g(q) = vel_g(q) + rotation_velocity(q) ! add the rotational velocity 
289289                            end do 
290290                        end if 
291291                    end if 
0 commit comments