import math
#from python_matrix_operations import *
import ctypes
import os
import sys
import numpy as np
# import elasticity_toolbox as elast
class ParameterSet(dict):
def __init__(self, *args, **kwargs):
super(ParameterSet, self).__init__(*args, **kwargs)
self.__dict__ = self
parameterSet = ParameterSet()
#---------------------------------------------------------------
#############################################
# Paths
#############################################
# Path for results and logfile
parameterSet.outputPath='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/perforated-bilayer/results'
parameterSet.baseName= 'perforated_wood_lower' #(needed for Output-Filename)
# Path for material description
# parameterSet.geometryFunctionPath =experiment/wood-bilayer/
#---------------------------------------------------------------
# Wooden bilayer, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
#--- define indicator function
# x[0] : y1-component -1/2 to 1/2
# x[1] : y2-component -1/2 to 1/2
# x[2] : x3-component range -1/2 to 1/2
#--- define indicator function
def indicatorFunction(x):
pRadius = math.sqrt(param_beta/np.pi)
if (x[2]>=(0.5-param_r)):
# if(np.sqrt(x[0]**2 + x[1]**2) < pRadius): #inside perforation
return 1 #Phase1
else :
if((x[0]**2 + x[1]**2) < pRadius**2): #inside perforation
return 3 #Phase3
else:
return 2 #Phase2
# def indicatorFunction(x):
# factor=1
# pRadius = 0.25
# if (x[2]>=(0.5-param_r) and np.sqrt(x[0]**2 + x[1]**2) < pRadius):
# return 3
# elif((x[2]>=(0.5-param_r))):
# return 1 #Phase1
# else :
# return 2 #Phase2
# # --- Number of material phases
# parameterSet.Phases=3
# def indicatorFunction(x):
# factor=1
# pRadius = 1
# # if (np.sqrt(x[0]*x[0] + x[1]*x[1]) < pRadius):
# if ((x[0] < 0 and math.sqrt(pow(x[1],2) + pow(x[2],2) ) < pRadius/4.0) or ( 0 < x[0] and math.sqrt(pow(x[1],2) + pow(x[2],2) ) > pRadius/4.0)):
# return 1
# else :
# return 2 #Phase2
# --- Number of material phases
parameterSet.Phases=3
# Parameters of the model
# -- (thickness upper layer) / (thickness)
# param_r = 0.22
param_r = 0.49
# -- thickness [meter]
param_h = 0.008
# -- moisture content in the flat state [%]
param_omega_flat = 17.01520754
# -- moisture content in the target state [%]
param_omega_target = 9.341730605
# -- Drehwinkel
param_theta = 0.0
# Design Parameter ratio between perforaton (cylindrical) volume and volume of upper layer
param_beta = 0.3
#
#
#
# -- increment of the moisture content
delta_omega=param_omega_target-param_omega_flat
# moisture content for material law
omega=param_omega_target
# --- Material properties from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6191116/#pone.0205607.ref015
# --- for European beech, moisture content omega = 15%
# --- L=direction orthogonal to layering and fibres = orthogonal to wood stem cross-section
# --- T=tangential zu layering
# --- R=orthogonal zu layering
# --- in MPa
# --- Properties are defined by affine function in dependence of moisture content omega via property = b_0+b_1 \omega
# --- coefficients of affine function are contained in the following array
# --- data taken from http://dx.doi.org/10.1016/j.cma.2014.10.031
properties_coefficients=np.array([
# [b_0, b_1]
[2565.6,-59.7], # E_R [MPa]
[885.4, -23.4], # E_T [MPa]
[17136.7,-282.4], # E_L [MPa]
[667.8, -15.19], # G_RT [MPa]
[1482, -15.26], # G_RL [MPa]
[1100, -17.72], # G_TL [MPa]
[0.2933, -0.001012], # nu_TR [1]
[0.383, -0.008722], # nu_LR [1]
[0.3368, -0.009071] # nu_LT [1]
])
# Compute actual material properties
E_R = properties_coefficients[0,0]+properties_coefficients[0,1]*omega
E_T = properties_coefficients[1,0]+properties_coefficients[1,1]*omega
E_L = properties_coefficients[2,0]+properties_coefficients[2,1]*omega
G_RT = properties_coefficients[3,0]+properties_coefficients[3,1]*omega
G_LR = properties_coefficients[4,0]+properties_coefficients[4,1]*omega
G_LT = properties_coefficients[5,0]+properties_coefficients[5,1]*omega
nu_TR = properties_coefficients[6,0]+properties_coefficients[6,1]*omega
nu_LR = properties_coefficients[7,0]+properties_coefficients[7,1]*omega
nu_LT = properties_coefficients[8,0]+properties_coefficients[8,1]*omega
# Compute the remaining Poisson ratios
nu_TL=nu_LT*E_T/E_L
nu_RT=nu_TR*E_R/E_T
nu_RL=nu_LR*E_R/E_L
#
# --- differential swelling strain
# --- relation to swelling strain eps: eps=alpha* delta_omega with delta_omega = change of water content in %
alpha_L=0.00011 # PLOS paper
alpha_R=0.00191 # PLOS paper
alpha_T=0.00462 # PLOS paper
# Umrechnen
#alpha_L=(1-1/(1+delta_omega*alpha_L))/delta_omega
#alpha_R=(1-1/(1+delta_omega*alpha_R))/delta_omega
#alpha_T=(1-1/(1+delta_omega*alpha_T))/delta_omega
# --- define geometry
# --- PHASE 1
# y_1-direction: L
# y_2-direction: T
# x_3-direction: R
# phase1_type="orthotropic"
# materialParameters_phase1 = [E_L,E_T,E_R,G_TL,G_RT,G_RL,nu_LT,nu_LR,nu_TR]
parameterSet.phase1_type="general_anisotropic"
[E_1,E_2,E_3]=[E_L,E_T,E_R]
[nu_12,nu_13,nu_23]=[nu_LT,nu_LR,nu_TR]
[nu_21,nu_31,nu_32]=[nu_TL,nu_RL,nu_RT]
[G_12,G_31,G_23]=[G_LT,G_LR,G_RT]
compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
[-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
[-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
materialParameters_phase1 = compliance_S
def prestrain_phase1(x):
# hB=delta_omega * alpha with delta_omega increment of moisture content and alpha swelling factor.
return [[1/param_h*delta_omega*alpha_L, 0, 0], [0,1/param_h*delta_omega*alpha_T,0], [0,0,1/param_h*delta_omega*alpha_R]]
# --- PHASE 2
# y_1-direction: R
# y_2-direction: L
# x_3-direction: T
parameterSet.phase2_type="general_anisotropic"
[E_1,E_2,E_3]=[E_R,E_L,E_T]
[nu_12,nu_13,nu_23]=[nu_RL,nu_RT,nu_LT]
[nu_21,nu_31,nu_32]=[nu_LR,nu_TR,nu_TL]
[G_12,G_31,G_23]=[G_LR,G_RT,G_LT]
compliance_S=np.array([[1/E_1, -nu_21/E_2, -nu_31/E_3, 0.0, 0.0, 0.0],
[-nu_12/E_1, 1/E_2, -nu_32/E_3, 0.0, 0.0, 0.0],
[-nu_13/E_1, -nu_23/E_2, 1/E_3, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1/G_23, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1/G_31, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 1/G_12]]);
materialParameters_phase2 = compliance_S
def prestrain_phase2(x):
return [[1/param_h*delta_omega*alpha_R, 0, 0], [0,1/param_h*delta_omega*alpha_L,0], [0,0,1/param_h*delta_omega*alpha_T]]
#Rotation um 2. Achse (= L)
parameterSet.phase2_axis = 1
# phase2_angle = param_theta
# -- Drehwinkel
parameterSet.phase2_angle = param_theta
# --- PHASE 3
parameterSet.phase3_type="isotropic"
epsilon = 1e-8
materialParameters_phase3 = [epsilon, epsilon]
def prestrain_phase3(x):
return [[0, 0, 0], [0,0,0], [0,0,0]]
# # --- PHASE 3 = Phase 1 gedreht
# # y_1-direction: L
# # y_2-direction: R
# # x_3-direction: T
# parameterSet.phase3_type="general_anisotropic"
# # Drehung um theta um Achse 2 = x_3-Achse
# N=elast.rotation_matrix_compliance(2,param_theta)
# materialParameters_phase3 = np.dot(np.dot(N,materialParameters_phase1),N.T)
# materialParameters_phase3 = 0.5*(materialParameters_phase3.T+materialParameters_phase3)
# # rotation of strain
# def prestrain_phase3(x):
# return elast.voigt_to_strain(np.dot(elast.rotation_matrix_compliance(2,param_theta),np.dot(elast.strain_to_voigt(np.array(prestrain_phase1(x))),N.T))).tolist()
# --- Choose scale ratio gamma:
parameterSet.gamma=1.0
#############################################
# Grid parameters
#############################################
## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
parameterSet.numLevels= '3 3' # computes all levels from first to second entry
# parameterSet.numLevels= '4 4'
#############################################
# Assembly options
#############################################
parameterSet.set_IntegralZero = 1 #(default = false)
parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
#parameterSet.arbitraryElementNumber = 3 #(default = 0)
#############################################
# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
#############################################
parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
#############################################
# Write/Output options #(default=false)
#############################################
# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
# --- (Additional debug output)
parameterSet.print_debug = 0 #(default=false)
# --- Write Correctos to VTK-File:
parameterSet.write_VTK = 1
# The grid can be refined several times for a higher resolution in the VTK-file.
parameterSet.subsamplingRefinement = 2
# --- (Optional output) L2Error, integral mean:
#parameterSet.write_L2Error = 1
#parameterSet.write_IntegralMean = 1
# --- check orthogonality (75) from paper:
parameterSet.write_checkOrthogonality = 0
# --- Write corrector-coefficients to log-File:
#parameterSet.write_corrector_phi1 = 1
#parameterSet.write_corrector_phi2 = 1
#parameterSet.write_corrector_phi3 = 1
# --- Print Condition number of matrix (can be expensive):
#parameterSet.print_conditionNumber= 1 #(default=false)
# --- write effective quantities to Matlab-folder for symbolic minimization:
parameterSet.write_toMATLAB = 1 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt