PAC_caseI.py

import subprocess
import re
import os
import numpy as np
import matplotlib.pyplot as plt
import math
import fileinput
import time
import matplotlib.ticker as tickers
import matplotlib as mpl
from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
import codecs
import sys
import threading

# Schreibe input datei für Parameter
def SetParametersCellProblem(ParameterSet, ParsetFilePath, outputPath):
    print('----set Parameters -----')
    with open(ParsetFilePath, 'r') as file:
        filedata = file.read()
        filedata = re.sub('(?m)^materialFunction\s?=.*','materialFunction = '+str(ParameterSet.materialFunction),filedata)
        filedata = re.sub('(?m)^gamma\s?=.*','gamma='+str(ParameterSet.gamma),filedata)
        filedata = re.sub('(?m)^numLevels\s?=.*','numLevels='+str(ParameterSet.numLevels)+' '+str(ParameterSet.numLevels) ,filedata)
        filedata = re.sub('(?m)^outputPath\s?=\s?.*','outputPath='+str(outputPath),filedata)
        f = open(ParsetFilePath,'w')
        f.write(filedata)
        f.close()


# Ändere Parameter der MaterialFunction
def SetParameterMaterialFunction(materialFunction, parameterName, parameterValue):
    with open(Path+"/"+materialFunction+'.py', 'r') as file:
        filedata = file.read()
        filedata = re.sub('(?m)^'+str(parameterName)+'\s?=.*',str(parameterName)+' = '+str(parameterValue),filedata)
        f = open(Path+"/"+materialFunction+'.py','w')
        f.write(filedata)
        f.close()

# Rufe Programm zum Lösen des Cell-Problems auf
def run_CellProblem(executable, parset,write_LOG):
    print('----- RUN Cell-Problem ----')
    processList = []
    LOGFILE = "Cell-Problem_output.log"
    print('LOGFILE:',LOGFILE)
    print('executable:',executable)
    if write_LOG:
        p = subprocess.Popen(executable + parset
                                        + " | tee " + LOGFILE, shell=True)

    else:
        p = subprocess.Popen(executable + parset
                                        + " | tee " + LOGFILE, shell=True)

    p.wait() # wait
    processList.append(p)
    exit_codes = [p.wait() for p in processList]

    return

# Read effective quantites
def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
    # Read Output Matrices (effective quantities)
    # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
    # -- Read Matrix Qhom
    X = []
    # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
    with codecs.open(QFilePath, encoding='utf-8-sig') as f:
        for line in f:
            s = line.split()
            X.append([float(s[i]) for i in range(len(s))])
    Q = np.array([[X[0][2], X[1][2], X[2][2]],
                  [X[3][2], X[4][2], X[5][2]],
                  [X[6][2], X[7][2], X[8][2]] ])

    # -- Read Beff (as Vector)
    X = []
    # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
    with codecs.open(BFilePath, encoding='utf-8-sig') as f:
        for line in f:
            s = line.split()
            X.append([float(s[i]) for i in range(len(s))])
    B = np.array([X[0][2], X[1][2], X[2][2]])
    return Q, B

# Function for evaluating the energy in terms of kappa, alpha and Q, B
def eval_energy(kappa,alpha,Q,B)  :
    G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
    return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]

#-------------------------------------------------------------------------------------------------------
########################
#### SET PARAMETERS ####
########################
# ----- Setup Paths -----
Path = "./experiment/PAC"
# parset = ' ./experiment/wood-bilayer/cellsolver.parset.wood'
ParsetFile = Path + '/cellsolver.parset'
executable = 'build-cmake/src/Cell-Problem'
write_LOG = True   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"

# ---------------------------------
# Setup Experiment
# ---------------------------------
outputPath = Path + '/results_caseI/'

# ----- Define Input parameters  --------------------
class ParameterSet:
    pass

ParameterSet.materialFunction  = "PAC"
ParameterSet.numLevels=4
case_=1
if case_==1:
    ParameterSet.vol = 0.16 # nu_f
    ParameterSet.t = 0.6 # R radius of fibre
    ParameterSet.r = 0.125 # R radius of fibre
    ParameterSet.h = 0.3
elif case_==2:
    ParameterSet.vol = 0.10 # nu_f
    ParameterSet.r = 0.1 # R radius of fibre
    ParameterSet.h = 0.3
    ParameterSet.t = 0.6 # t thickness of hinge,
elif case_==3:
    ParameterSet.vol = 0.16 # nu_f
    ParameterSet.r = 0.1 # R radius of fibre
    ParameterSet.h = 0.2
    ParameterSet.t = 0.6 # t thickness of hinge,
# 
vol_=np.pi*ParameterSet.r**2
# h * epsilon_ * vol = vol_
epsilon_ = vol_/(ParameterSet.h*ParameterSet.vol)
# gamma
ParameterSet.gamma = ParameterSet.t/epsilon_

# ----- Define Parameters for Material Function  --------------------
# Liste mit Drehwinkel eigenstraint
materialFunctionParameter=[0.1, 0.2, 0.3]

# ------ Loops through Parameters for Material Function -----------
for i in range(0,np.shape(materialFunctionParameter)[0]):
    print("------------------")
    print("New Loop")
    print("------------------")
   # Check output directory
    path = outputPath + str(i)
    isExist = os.path.exists(path)
    if not isExist:
        # Create a new directory because it does not existl
        os.makedirs(path)
        print("The new directory " + path + " is created!")
    # keine Parameter daher naechste Zeiel auskommentiert
    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_eigenstrain",materialFunctionParameter[i])    
    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_t", ParameterSet.t)    
    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_h", ParameterSet.h)    
    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_vol", ParameterSet.vol)    
    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_r", ParameterSet.r)    
    SetParameterMaterialFunction(ParameterSet.materialFunction, "param_gamma", ParameterSet.gamma)    
    SetParametersCellProblem(ParameterSet, ParsetFile, path)
    #Run Cell-Problem
    thread = threading.Thread(target=run_CellProblem(executable, " ./"+ParsetFile, write_LOG))
    thread.start()
    # ---------------------------------------------------
    # wait here for the result to be available before continuing
    thread.join()
    f = open(path+"/parameter.txt", "w")
    f.write("param_eigenstrain = "+str(materialFunctionParameter[i])+"\n")
    f.close()   
    #