-
Jan Habscheid authoredJan Habscheid authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
GenerateMicrostructures_fun.py 3.35 KiB
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