Skip to content
Snippets Groups Projects
perforated_wood_upper.py 10.59 KiB
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_upper'   #(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 
        if((x[0]**2 + x[1]**2) < pRadius**2):  #inside perforation     
            return 3  #Phase3
        else:  
            return 1  #Phase1
    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