Add Model Patches

Author: henryleberre

docs/documentation/case.md

@@ -0,0 +1,84 @@
+import math, argparse
+
+from PIL import Image, ImageFilter
+
+parser = argparse.ArgumentParser(prog='img2stl', description='Turn an image into an STL file')
+parser.add_argument('input',  type=str, help='Path to the image file')
+parser.add_argument('output', type=str, help='Path to the output OBJ file')
+
+VARS = vars(parser.parse_args())
+
+img = Image.open(VARS['input'])
+
+W, H = img.size
+
+pixels = img.load()
+for y in range(H):
+    for x in range(W):
+        if pixels[x,y][3] < 255:
+            pixels[x,y] = (0,0,0,255)
+
+img.save('temp.png')
+img = img.convert('L')
+
+THRESHOLD = 0.1
+
+with open(VARS['output'], 'w') as f:
+    i = 0
+    pixels = img.load()
+    for y in range(H):
+        for x in range(W):
+            c = pixels[x,y]
+
+            if c < 255 * THRESHOLD: continue
+
+            area = c / 255.0
+            d    = (1.0 - math.sqrt(area)) / 2.0
+
+            def tx(_x: int, _f: bool):
+                if _f: _x = _x + d*(1 - _x*2)
+                return -1.0+2.0 * (x + _x)/float(W)
+
+            def ty(_y: int, _f: bool):
+                if _f: _y = _y + d*(1 - _y*2)
+                return  1.0-2.0 * ((1-H/float(W))/2 + (y+_y)/float(W))
+
+            def tz(_z: int, _f: bool):
+                if _f: _z = _z + d*(1 - _z*2)
+                return max(1.0 / float(W), 1.0 / float(H)) * (-1.0+2.0 * _z)
+
+            def p1(_x: int, _y: int, _z: int):
+                return tx(_x, True), ty(_y, True), tz(_z, False)
+
+            for _x, _y, _z in [p1(0,0,0), p1(1,0,0), p1(0,1,0), p1(1,1,0),
+                               p1(0,0,1), p1(1,0,1), p1(0,1,1), p1(1,1,1)]:
+                f.write(f'v {_x} {_y} {_z}\n')
+            
+            f.write(f'f {i+1} {i+2} {i+3}\n'); f.write(f'f {i+2} {i+3} {i+4}\n')
+            f.write(f'f {i+5} {i+6} {i+7}\n'); f.write(f'f {i+6} {i+7} {i+8}\n')
+
+            i = i + 8
+
+            def p2(_x: int, _y: int, _z: int):
+                return tx(_x, False), ty(_y, True), tz(_z, True)
+
+            for _x, _y, _z in [p2(0,0,0), p2(1,0,0), p2(0,1,0), p2(1,1,0),
+                               p2(0,0,1), p2(1,0,1), p2(0,1,1), p2(1,1,1)]:
+                f.write(f'v {_x} {_y} {_z}\n')
+
+            f.write(f'f {i+1} {i+5} {i+7}\n'); f.write(f'f {i+1} {i+3} {i+7}\n')
+            f.write(f'f {i+2} {i+6} {i+8}\n'); f.write(f'f {i+2} {i+4} {i+8}\n')
+
+            i = i + 8
+
+            def p3(_x: int, _y: int, _z: int):
+                return tx(_x, True), ty(_y, False), tz(_z, True)
+
+            for _x, _y, _z in [p3(0,0,0), p3(1,0,0), p3(0,1,0), p3(1,1,0),
+                               p3(0,0,1), p3(1,0,1), p3(0,1,1), p3(1,1,1)]:
+                f.write(f'v {_x} {_y} {_z}\n')
+
+            f.write(f'f {i+1} {i+2} {i+5}\n'); f.write(f'f {i+2} {i+5} {i+6}\n')
+            f.write(f'f {i+3} {i+4} {i+7}\n'); f.write(f'f {i+4} {i+7} {i+8}\n')
+
+            i = i + 8

misc/img2stl.py

@@ -2,6 +2,8 @@
 !! @file m_derived_types.f90
 !! @brief Contains module m_derived_types
 
+#:include "macros.fpp"
+
 !> @brief This file contains the definitions of all of the custom-defined
 !!              types used in the pre-process code.
 module m_derived_types
@@ -59,6 +61,45 @@ module m_derived_types
         integer, dimension(:, :, :), allocatable :: fullmom !< Moment indices for qbmm
     end type bub_bounds_info
 
+    !> Defines parameters for a Model Patch
+    type :: ic_model_parameters
+        character(LEN=pathlen_max) :: filepath !<
+        !! Path the STL file relative to case_dir.
+
+        t_vec3 :: translate !<
+        !! Translation of the STL object.
+
+        t_vec3 :: scale !<
+        !! Scale factor for the STL object.
+
+        t_vec3 :: rotate !<
+        !! Angle to rotate the STL object along each cartesian coordinate axis,
+        !! in radians.
+
+        integer :: spc !<
+        !! Number of samples per cell to use when discretizing the STL object.
+    end type ic_model_parameters
+
+    type :: t_triangle
+        real(kind(0d0)), dimension(1:3, 1:3) :: v ! Vertices of the triangle
+        t_vec3 :: n ! Normal vector
+    end type t_triangle
+
+    type :: t_ray
+        t_vec3 :: o ! Origin
+        t_vec3 :: d ! Direction
+    end type t_ray
+
+    type :: t_bbox
+        t_vec3 :: min ! Minimum coordinates
+        t_vec3 :: max ! Maximum coordinates
+    end type t_bbox
+
+    type :: t_model
+        integer                       :: ntrs   ! Number of triangles
+        type(t_triangle), allocatable :: trs(:) ! Triangles
+    end type t_model
+
     !> Derived type adding initial condition (ic) patch parameters as attributes
     !! NOTE: The requirements for the specification of the above parameters
     !! are strongly dependent on both the choice of the multicomponent flow
@@ -79,6 +120,8 @@ module m_derived_types
         !! patch geometries. It is specified through its x-, y-, and z-components
         !! respectively.
 
+        type(ic_model_parameters) :: model !< Model parameters
+
         real(kind(0d0)) :: epsilon, beta !<
         !! The isentropic vortex parameters administrating, respectively, both
         !! the amplitude of the disturbance as well as its domain of influence.

src/common/m_derived_types.f90 renamed to src/common/m_derived_types.fpp

@@ -2,6 +2,7 @@
 !>
 !! @file m_helper.f90
 !! @brief Contains module m_helper
+
 module m_helper
 
     ! Dependencies =============================================================
@@ -20,7 +21,14 @@ module m_helper
         s_initialize_nonpoly, &
         s_simpson, &
         s_transcoeff, &
-        s_int_to_str
+        s_int_to_str, &
+        s_transform_vec, &
+        s_transform_triangle, &
+        s_transform_model, &
+        s_swap, &
+        f_cross, &
+        f_create_transform_matrix, &
+        f_create_bbox
 
 contains
 
@@ -148,15 +156,9 @@ contains
         rhol0 = rhoref
         pl0 = pref
 
-#ifdef MFC_SIMULATION
         @:ALLOCATE(pb0(nb), mass_n0(nb), mass_v0(nb), Pe_T(nb))
         @:ALLOCATE(k_n(nb), k_v(nb), omegaN(nb))
         @:ALLOCATE(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))
-#else
-        allocate(pb0(nb), mass_n0(nb), mass_v0(nb), Pe_T(nb))
-        allocate(k_n(nb), k_v(nb), omegaN(nb))
-        allocate(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))
-#endif
 
         pb0(:) = dflt_real
         mass_n0(:) = dflt_real
@@ -331,5 +333,150 @@ contains
         tmp = DEXP(-0.5d0*(phi(nb)/sd)**2)/DSQRT(2.d0*pi)/sd
         weight(nb) = tmp*dphi/3.d0
     end subroutine s_simpson
+    
+    !> This procedure computes the cross product of two vectors.
+    !! @param a First vector.
+    !! @param b Second vector.
+    !! @return The cross product of the two vectors.
+    function f_cross(a, b) result(c)
+        real(kind(0d0)), dimension(3), intent(in) :: a, b
+        real(kind(0d0)), dimension(3)             :: c
+
+        c(1) = a(2) * b(3) - a(3) * b(2)
+        c(2) = a(3) * b(1) - a(1) * b(3)
+        c(3) = a(1) * b(2) - a(2) * b(1)
+    end function f_cross
+
+    !> This procedure swaps two real numbers.
+    !! @param lhs Left-hand side.
+    !! @param rhs Right-hand side.
+    subroutine s_swap(lhs, rhs)
+        real(kind(0d0)), intent(inout) :: lhs, rhs
+        real(kind(0d0))                :: ltemp
+
+        ltemp = lhs
+        lhs   = rhs
+        rhs   = ltemp
+    end subroutine s_swap
+
+    !> This procedure creates a transformation matrix.
+    !! @param  p Parameters for the transformation.
+    !! @return Transformation matrix.
+    function f_create_transform_matrix(p) result(out_matrix)
+
+        type(ic_model_parameters) :: p
+
+        t_mat4x4 :: sc, rz, rx, ry, tr, out_matrix
+
+        sc = transpose(reshape([ &
+            p%scale(1), 0d0,        0d0,        0d0, &
+            0d0,        p%scale(2), 0d0,        0d0, &
+            0d0,        0d0,        p%scale(3), 0d0, &
+            0d0,        0d0,        0d0,        1d0  ], shape(sc)))
+
+        rz = transpose(reshape([ &
+            cos(p%rotate(3)), -sin(p%rotate(3)), 0d0, 0d0, &
+            sin(p%rotate(3)),  cos(p%rotate(3)), 0d0, 0d0, &
+            0d0,               0d0,              1d0, 0d0, &
+            0d0,               0d0,              0d0, 1d0  ], shape(rz)))
+
+        rx = transpose(reshape([ &
+            1d0, 0d0,               0d0,              0d0, &
+            0d0, cos(p%rotate(1)), -sin(p%rotate(1)), 0d0, &
+            0d0, sin(p%rotate(1)),  cos(p%rotate(1)), 0d0, &
+            0d0, 0d0,               0d0,              1d0  ], shape(rx)))
+
+        ry = transpose(reshape([ &
+            cos(p%rotate(2)),  0d0, sin(p%rotate(2)), 0d0, &
+            0d0,               1d0, 0d0,              0d0, &
+            -sin(p%rotate(2)), 0d0, cos(p%rotate(2)), 0d0, &
+            0d0,               0d0, 0d0,              1d0  ], shape(ry)))
+
+        tr = transpose(reshape([ &
+            1d0, 0d0, 0d0, p%translate(1), &
+            0d0, 1d0, 0d0, p%translate(2), &
+            0d0, 0d0, 1d0, p%translate(3), &
+            0d0, 0d0, 0d0, 1d0 ], shape(tr)))
+
+        out_matrix = matmul(tr, matmul(ry, matmul(rx, matmul(rz, sc))))
+
+    end function f_create_transform_matrix
+
+    !> This procedure transforms a vector by a matrix.
+    !! @param vec Vector to transform.
+    !! @param matrix Transformation matrix.
+    subroutine s_transform_vec(vec, matrix)
+
+        t_vec3,   intent(inout) :: vec
+        t_mat4x4, intent(in)    :: matrix
+
+        real(kind(0d0)), dimension(1:4) :: tmp
+
+        tmp = matmul(matrix, [ vec(1), vec(2), vec(3), 1d0 ])
+        vec = tmp(1:3)
+
+    end subroutine s_transform_vec
+
+    !> This procedure transforms a triangle by a matrix, one vertex at a time.
+    !! @param triangle Triangle to transform.
+    !! @param matrix   Transformation matrix.
+    subroutine s_transform_triangle(triangle, matrix)
+
+        type(t_triangle), intent(inout) :: triangle
+        t_mat4x4, intent(in) :: matrix
+
+        integer :: i
+
+        real(kind(0d0)), dimension(1:4) :: tmp
+
+        do i = 1, 3
+            call s_transform_vec(triangle%v(i,:), matrix)
+        end do
+
+    end subroutine s_transform_triangle
+
+    !> This procedure transforms a model by a matrix, one triangle at a time.
+    !! @param model  Model to transform.
+    !! @param matrix Transformation matrix.
+    subroutine s_transform_model(model, matrix)
+
+        type(t_model), intent(inout) :: model
+        t_mat4x4,      intent(in)    :: matrix
+
+        integer :: i
+
+        do i = 1, size(model%trs)
+            call s_transform_triangle(model%trs(i), matrix)
+        end do
+
+    end subroutine s_transform_model
+
+    !> This procedure creates a bounding box for a model.
+    !! @param model Model to create bounding box for.
+    !! @return Bounding box.
+    function f_create_bbox(model) result(bbox)
+
+        type(t_model), intent(in) :: model
+        type(t_bbox)              :: bbox
+
+        integer :: i, j
+
+        if (size(model%trs) == 0) then
+            bbox%min = 0d0
+            bbox%max = 0d0
+            return
+        end if
+
+        bbox%min = model%trs(1)%v(1,:)
+        bbox%max = model%trs(1)%v(1,:)
+
+        do i = 1, size(model%trs)
+            do j = 1, 3
+                bbox%min = min(bbox%min, model%trs(i)%v(j,:))
+                bbox%max = max(bbox%max, model%trs(i)%v(j,:))
+            end do
+        end do
+
+    end function f_create_bbox
 
 end module m_helper

src/common/m_helper.fpp

@@ -132,7 +132,7 @@ contains
         ! Depending on model_eqns and bubbles, the appropriate procedure
         ! for computing pressure is targeted by the procedure pointer
 
-        if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then   
+        if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then  
             pres = (energy - dyn_p - pi_inf)/gamma
         else if ((model_eqns /= 4) .and. bubbles) then
             pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf)/gamma

src/common/m_variables_conversion.fpp

@@ -79,7 +79,9 @@ contains
                 elseif (patch_icpp(i)%geometry == 19) then
                     print *, '3d var circle'
                 elseif (patch_icpp(i)%geometry == 20) then
-                    call s_check_2D_TaylorGreen_vortex_patch_geometry(i)    
+                    call s_check_2D_TaylorGreen_vortex_patch_geometry(i)
+                elseif (patch_icpp(i)%geometry == 21) then
+                    call s_check_model_geometry(i)
                 else
                     call s_mpi_abort('Unsupported choice of the '// &
                             'geometry of active patch '//trim(iStr)//&
@@ -851,4 +853,23 @@ contains
 
     end subroutine s_check_inactive_patch_primitive_variables ! ------------
 
+    subroutine s_check_model_geometry(patch_id) ! ------------------------------
+
+        integer, intent(IN) :: patch_id
+
+        logical :: file_exists
+
+        inquire(file=patch_icpp(patch_id)%model%filepath, exist=file_exists)
+
+        if (.not. file_exists) then
+
+            print '(A,I0,A)', 'Model file '//trim(patch_icpp(patch_id)%model%filepath)//&
+                ' requested by patch ', patch_id,' does not exist. Exiting ...'
+
+            call s_mpi_abort()
+
+        end if
+
+    end subroutine s_check_model_geometry ! -----------------------------------
+
 end module m_check_patches

src/pre_process/m_check_patches.fpp

@@ -26,6 +26,14 @@ subroutine s_check_inputs()
         logical :: dir_check !<
             !! Logical variable used to test the existence of folders
 
+#ifndef MFC_MPI
+        if (parallel_io .eqv. .true.) then
+            print '(A)', 'MFC built with --no-mpi requires parallel_io=F. ' // &
+                'Exiting ...'
+            call s_mpi_abort()
+        end if
+#endif
+
         bub_fac = 0
         if (bubbles .and. (num_fluids == 1)) bub_fac = 1
         ! Startup checks for bubbles and bubble variables