GenerateMicrostructures_fun.py

import numpy as np
# import scipy.optimize as scopt
import splinepy as sp
import pandas as pd

from gismo.gismo_export import AdditionalBlocks, export

def generate(size_1, size_2, size_3, BOX_LENGTH, BOX_HEIGHT, EPS, INLET_BOUNDARY_ID, OUTLET_BOUNDARY_ID, knots_y, TILING, CLOSING_FACE, MICROTILE):
    # Define microstructure deformation function
    initial_macro_spline = sp.helpme.create.box(BOX_LENGTH, BOX_HEIGHT)
    # x0_max, x1_max = initial_macro_spline.cps.max(axis=0)
    # print(f'x0_max: {x0_max}, x1_max: {x1_max}')
    
    # Define identifier functions for microstructure boundaries
    def identifier_inlet(points):
        return points[:, 0] < EPS

    def identifier_outlet(points):
        return points[:, 0] > BOX_LENGTH - EPS

    boundary_identifier_dict = {
        identifier_inlet: INLET_BOUNDARY_ID,
        identifier_outlet: OUTLET_BOUNDARY_ID,
    }

    # Define current parameterization function
    def para_function_(x):
        # print(f'x: {x}')
        # print(f'index: {index}')
        x_ = x[0][0]
        y_ = x[0][1]
        if np.isclose(y_, knots_y[0], rtol=1e-2):
            return np.array(size_1).reshape(-1, 1)
        elif np.isclose(y_, knots_y[1], rtol=1e-2):
            return np.array(size_2).reshape(-1, 1)
        elif np.isclose(y_, knots_y[2], rtol=1e-2):
            return np.array(size_3).reshape(-1, 1)
        else:
            raise ValueError(f'Invalid y_: {y_}')

    # Create microstructure
    generator = sp.microstructure.Microstructure(
        deformation_function=initial_macro_spline.copy(),
        # deformation_function=initial_parameter_spline.copy(),
        tiling=TILING,
        microtile=MICROTILE,
        parametrization_function=para_function_,
    )

    microstructure = generator.create(
        closing_face=CLOSING_FACE,
    )

    # Simulation parameters
    for (
        identifier_function,
        boundary_id,
    ) in boundary_identifier_dict.items():
        microstructure.boundary_from_function(
            identifier_function, boundary_id=boundary_id
        )

    # additional_blocks = sp.io.gismo.AdditionalBlocks()
    additional_blocks = AdditionalBlocks()
    # Velocity boundary conditions
    additional_blocks.add_boundary_conditions(
        block_id=2,
        dim=2,
        function_list=["0", (f"y * ({BOX_HEIGHT}-y)", "0")],
        bc_list=[
            (f"BID{INLET_BOUNDARY_ID}", "Dirichlet", 1),  # Inlet
            ("BID1", "Dirichlet", 0),  # Walls
        ],
        unknown_id=1,
        multipatch_id=0,
        comment=" Velocity boundary conditions: parabolic inflow field ",
    )

    # Pressure boundary conditions
    additional_blocks.add_boundary_conditions(
        block_id=3,
        dim=2,
        function_list=["0"],
        # bc_list=[(f"BID{OUTLET_BOUNDARY_ID}", "Dirichlet", 0)],
        bc_list=[], # use empty list, as bc_list argument is needed, but cv_list is wanted for corner value only
        cv_list=[
            ("0", "40", "2", "0"),
        ],  # unknown, patch, corner, function
        unknown_id=0,
        multipatch_id=0,
        comment=" Pressure boundary conditions: fix pressure to zero at the corner value",
    )

    # Get default assembly options
    additional_blocks.add_assembly_options(
        block_id=4, comment=" Assembler options "
    )

    gismo_export_options = additional_blocks.to_list()


    return microstructure, gismo_export_options