Skip to content
Snippets Groups Projects
Commit a4c0517f authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Add microstructure_testsuite

parent 78d8e89a
No related branches found
No related tags found
No related merge requests found
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 sympy as sym
# from prettytable import PrettyTable
# import latextable
# from texttable import Texttable
# from tabulate import tabulate
#-------------------------------------------------------------------------------------------------------
def run_CellProblem(executable, parset, gridLevel, gamma, mu1,lambda1, rho1, alpha, beta, theta, material_prestrain_imp, outputPath, write_materialFunctions, write_prestrainFunctions, write_VTK,write_L2Error ,write_IntegralMean, write_LOG ):
print('----- RUN Cell-Problem ----')
processList = []
LOGFILE = "Cell-Problem_output.log"
print('LOGFILE:',LOGFILE)
print('executable:',executable)
print('parset:',parset)
# start_time = time.time()
if write_LOG:
p = subprocess.Popen(executable + parset
# + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
+ " -gamma " + str(gamma)
+ " -mu1 " + str(mu1)
+ " -lambda1 " + str(lambda1)
+ " -alpha " + str(alpha)
+ " -beta " + str(beta)
+ " -theta " + str(theta)
+ " -material_prestrain_imp " + str(material_prestrain_imp)
+ " -write_materialFunctions " + str(write_materialFunctions)
+ " -outputPath " + str(outputPath)
+ " | tee " + LOGFILE, shell=True)
else:
p = subprocess.Popen(executable + parset
# + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
+ " -gamma " + str(gamma)
+ " -mu1 " + str(mu1)
+ " -lambda1 " + str(lambda1)
+ " -alpha " + str(alpha)
+ " -beta " + str(beta)
+ " -theta " + str(theta)
+ " -material_prestrain_imp " + str(material_prestrain_imp)
+ " -write_materialFunctions " + str(write_materialFunctions)
+ " -outputPath " + str(outputPath), stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, shell=True ) # surpress Logging
p.wait() # wait
# print("--- %s seconds ---" % (time.time() - start_time))
# print('------FINISHED PROGRAM on level:' + str(gridLevel))
processList.append(p)
###Wait for all simulation subprocesses before proceeding
exit_codes = [p.wait() for p in processList]
return
#-------------------------------------------------------------------------------------------------------
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
#-------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------
########################
#### SET PARAMETERS ####
########################
# ----- Setup Paths -----
outputPath = " ../outputs"
QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt'
BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'
print('---- Input parameters: -----')
mu1 = 10.0
rho1 = 1.0
alpha = 2.8
beta = 2.0
theta = 1.0/4.0
gamma = 0.75
lambda1= 10.0
gridLevel = 4
#############################################
# Choose preferred Geometry/Prestrain/Material
#############################################
#material_prestrain_imp= "analytical_Example"
material_prestrain_imp= "parametrized_Laminate"
#############################################
# Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER
#############################################
Solvertype = 1 #(default = 1)
Solver_verbosity = 0 #(default = 2) degree of information for solver output
set_IntegralZero = True
#############################################
# Output-Options
#############################################
write_materialFunctions = False
write_prestrainFunctions = False
write_VTK = False
# write_L2Error = True
# write_IntegralMean = True
write_L2Error = False
write_IntegralMean = False
write_LOG = False # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
print('mu1: ', mu1)
print('rho1: ', rho1)
print('alpha: ', alpha)
print('beta: ', beta)
print('theta: ', theta)
print('gamma:', gamma)
print('material_prestrain_imp: ', material_prestrain_imp)
print('gridLevel: ', gridLevel)
parset = ' ../inputs/cellsolver.parset'
executable = ' ../build-cmake/src/Cell-Problem'
print('---------------------------------------------------------')
###########################################################################################################
run_CellProblem(executable,
parset,
gridLevel,
gamma,
mu1,
lambda1,
rho1,
alpha,
beta,
theta,
material_prestrain_imp,
outputPath,
write_materialFunctions,
write_prestrainFunctions,
write_VTK,
write_L2Error,
write_IntegralMean,
write_LOG
)
print('Read effective quantities...')
Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
# Q, B = ReadEffectiveQuantities()
print('Q:', Q)
print('B:', B)
#
#
# for kappa in [4]:
#
# print('--- RUN COMPUTATION FOR KAPPA=',kappa)
#
# # --- map R2 into R2
# if domainDim == 2 and targetDim ==2:
# print('domainDim == 2 and targetDim == 2')
#
# if kappa == 0 or kappa == 1:
# initialIterate = "initialIterateR2_R2"
# else:
# initialIterate = "initialIterateR2_R2_deg"+str(kappa)
# # initialIterate = "initialIterateR2_R2"
# # initialIterate = "initialIterateR2_R2_deg2"
# # initialIterate = "initialIterateR2_R2_deg3"
# # initialIterate = "initialIterateR2_R2_deg4"
# # initialIterate = "initialIterateR2_R2_deg5"
# # initialIterate = "initialIterateR2_R2_deg8"
#
# executable = '../build-cmake/src/harmonicmaps_2d'
# parset = ' ../parsets/harmonicmaps_EOC_2d.parset'
#
# # --- map R2 into R3
# if domainDim == 2 and targetDim ==3:
# print('domainDim == 2 and targetDim == 3')
#
# # initialIterate = "initialIterateR2_R3"
# # initialIterate = "initialIterateR2_R3_smoothed"
# if kappa == 0 or kappa == 1:
# # initialIterate = "initialIterateR2_R3"
# initialIterate = "initialIterateR2_R3_smoothed"
# else:
# initialIterate = "initialIterateR2_R2_deg"+str(kappa)
#
# # initialIterate = "initialIterateR2_R3_deg2"
# # initialIterate = "initialIterateR2_R3_deg3"
# # initialIterate = "initialIterateR2_R3_deg4"
# # initialIterate = "initialIterateR2_R3_deg5"
# executable = '../build-cmake/src/harmonicmaps_2d'
# parset = ' ../parsets/harmonicmaps_EOC_2d.parset'
#
# # --- map R3 into R3
# if domainDim == 3 and targetDim ==3:
# print('domainDim == 3 and targetDim == 3')
#
# if kappa == 0 or kappa == 1:
# initialIterate = "initialIterateR3_R3"
# else:
# initialIterate = "initialIterateR3_R3_deg"+str(kappa)
#
#
# # initialIterate = "initialIterateR3_R3_pertubed"
# # initialIterate = "initialIterateR2_R2"
# # initialIterate = "initialIterateR3_R3_deg2"
# # initialIterate = "initialIterateR3_R3_deg3"
# # initialIterate = "initialIterateR3_R3_deg4"
# # initialIterate = "initialIterateR3_R3_deg5"
# # initialIterate = "initialIterateR3_R3_deg8"
# # initialIterate = "initialIterateR3_R3_deg12"
# # initialIterate = "initialIterateR3_R3_deg3_Pert"
# # initialIterate = "initialIterateR3_R3_deg3_newOrigin"
#
#
# executable = '../build-cmake/src/harmonicmaps_3d'
# parset = ' ../parsets/harmonicmaps_EOC_3d.parset'
#
#
#
#
#
#
#
# print(' InitialIterate used: ', initialIterate)
#
# path = os.getcwd()
# print("Path: ", path)
#
# for order in [1]:
# minLevel = 5 + 1 - order;
# maxLevel = 10 + 1 - order;
# minLevel = 2 + 1 - order;
# maxLevel = 8 + 1 - order;
#--- Setup LatexTable
# table_1 = Texttable()
# # table_1.maketitle = "Test"
# table_1.set_cols_align(["l", "r", "c"])
# table_1.set_cols_valign(["t", "m", "b"])
# table_1.add_rows([["Name", "Age", "Nickname"],
# ["Mr\nXavier\nHuon", 32, "Xav'"],
# ["Mr\nBaptiste\nClement", 1, "Baby"],
# ["Mme\nLouise\nBourgeau", 28, "Lou\n \nLoue"]])
#
# # table_1.add_rows("\multicolumn{3}{c}{ Projection-based-FE of order: 1, $\kappa = 2$}")
# print('-- Example 1: Basic --')
# print('Texttable Output:')
# print(table_1.draw())
# print('\nLatextable Output:')
# print(latextable.draw_latex(table_1, caption="An example table.", label="table:example_table"))
#
# rows = [['Rocket', 'Organisation', 'LEO Payload (Tonnes)', 'Maiden Flight'],
# ['Saturn V', 'NASA', '140', '1967'],
# ['Space Shuttle', 'NASA', '24.4', '1981'],
# ['Falcon 9 FT-Expended', 'SpaceX', '22.8', '2017'],
# ['Ariane 5 ECA', 'ESA', '21', '2002']]
#
# table = Texttable()
# table.set_cols_align(["c"] * 4)
# table.set_deco(Texttable.HEADER | Texttable.VLINES)
# table.add_rows(rows)
#
# print('Tabulate Table:')
# print(tabulate(rows, headers='firstrow'))
#
# # print('\nTexttable Table:')
# # print(table.draw())
#
# print('\nTabulate Latex:')
# print(tabulate(rows, headers='firstrow', tablefmt='latex'))
# print('\nTexttable Latex:')
# print(latextable.draw_latex(table, caption="A comparison of rocket features."))
#
# #--- Setup Output-Table:
# x = PrettyTable()
# x.title = interpolationMethod + "-FE" + ' of order ' + str(order) + " into R" + str(targetDim)
# # x.field_names = ["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$", "$\lnorm{u_h-u_{2h}$", "$\lnorm{u_{2h}-u_{4h}$", "$\seminorm{u_h-u_{2h}}$", "$\seminorm{u_{2h}-u_{4h}}$" ]
# x.field_names = ["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$", "$\lnorm{u_h-u_{2h}$", "$\lnorm{u_{2h}-u_{4h}$", "$\seminorm{u_h-u_{2h}}$", "$\seminorm{u_{2h}-u_{4h}}$" , "wall-time" ]
# x.align["k"] = "l" # Left align city names
# x.padding_width = 1 # One space between column edges and contents (default)
#
# rows = []
# rows.append(["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$" , "wall-time" ])
# # rows.append(["k", "#Triang/#DOF", "# TR-Steps","wall-time","$E^{Di}(0)$", "$E^{Di}$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$" ])
# print('rows:', rows)
#
# if compute_HarmonicMaps:
#
# computeHarmonicMap(minLevel,maxLevel,domainDim,targetDim,order,interpolationMethod,maxTrustRegionSteps,initialIterate,randomInitialIterate,readConfiguration,configurationFile,pertubation,pertubationRadius,epsilon,tolerance,executable,parset )
#
# # --------------------
#
#
# if compute_Error:
# print("-------------------------------------------------------")
# print("Now measuring errors with discretizationErrorMode:" + str(discretizationErrorMode))
# # subprocess.call(["echo", "Now measuring errors"])
#
# # EOC_l2_list = ["-","-"]
#
# energyList = []
# initial_energyList = []
# stepList = []
# ndofList = []
# timeList = []
# # EOC_l2_list = []
#
# #
# for numLevels in range(minLevel,maxLevel+1):
#
# print("NUMBER OF LEVELS:", numLevels)
# # Measure the discretization errors against the solution on the finest grid
# # LOGFILE = "./compute-disc-error_" + str(order) + "_" + str(numLevels) + ".log"
#
#
# #subprocess.Popen(["../build-cmake/src/compute-disc-error", "compute-disc-error-skyrmions-hexagon.parset",
# #"-order", str(order),
# #"-numLevels", str(numLevels),
# #"-numReferenceLevels", str(maxLevel),
# #"-simulationData", "harmonicmaps-result-" + str(order) + "-" + str(numLevels) + ".data",
# #"-referenceData", "harmonicmaps-result-" + str(order) + "-" + str(maxLevel) + ".data"])
#
#
# # ------------ Get Energy / nDofs/Steps etc. -----------------------------------------------------------------------------
# LOGFILE_comp = "./harmonicmapsR"+ str(domainDim) +"_intoR"+ str(targetDim) + "_deg" + str(kappa) + "_"+ interpolationMethod + "_" + str(order) + "_" + str(numLevels) + ".log"
# # Read Energy Values:
# with open(LOGFILE_comp, 'r') as file:
# output = file.read()
#
# try:
# # tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1) # muss nun nichtmehr am Zeilenanfang stehen! :)
# # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1)
# # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# tmp_energy = re.findall(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
# tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
# tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
# # tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# # tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
# tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+\d?)',output)
#
# except AttributeError:
# # tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output) # muss nun nichtmehr am Zeilenanfang stehen! :)
# tmp_energy = re.search(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
# tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
# tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
# # tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# # tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
# tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+\d?)',output)
#
# if maxTrustRegionSteps>0 :
# print('tmp_step:',tmp_step)
# print('tmp_step last:', int(tmp_step[-1]))
#
# print('tmp_dofs:',tmp_dofs)
# print('tmp_dofs last:', int(tmp_dofs[-1]))
#
# print('type tmp_energy:', type(tmp_energy))
# print('tmp_energy:', tmp_energy)
# print('tmp_energy last:', tmp_energy[-1])
#
# print('tmp_time:', tmp_time)
#
# energy = float(tmp_energy[-1]) #[1] since otherwise it recognizes "2" from L2error...
# energyList.append(energy)
#
# initialEnergy = float(tmp_energy[0])
# initial_energyList.append(initialEnergy)
#
# steps = int(tmp_step[-1])+1 #starts with zero therefore +1
# stepList.append(steps)
#
# ndofs = int(tmp_dofs[-1])
# ndofList.append(ndofs)
#
# time = float(tmp_time[-1])
# timeList.append(time)
# # -----------------------------------------------------------------------------------------------
#
# if numLevels >= (minLevel+2) and discretizationErrorMode=='discrete':
#
# [EOC_l2, EOC_h1, EOC_energy, L2_fine , L2_coarse , H1_fine , H1_coarse, constraintError] \
# = computeError_discrete(numLevels,minLevel,domainDim,targetDim,order,interpolationMethod,discretizationErrorMode,maxTrustRegionSteps,energyList,ndofList,stepList,initial_energyList)
#
# # x.add_row([numLevels, currentDofs , currentSteps,currentInitialEnergy, currentEnergy, EOC_l2, EOC_h1, EOC_energy , format(l2_fine, "5.3e"), format(l2_coarse, "5.3e"), format(h1_fine, "5.3e"), format(h1_coarse, "5.3e") ])
# # x.add_row([numLevels, "-" , "-", "energy", EOC_l2, EOC_h1, "EOC_h^{E}" , format(l2_fine, "5.3e"), format(l2_coarse, "5.3e"), format(h1_fine, "5.3e"), format(h1_coarse, "5.3e") ])
#
# elif numLevels >= (minLevel+1) and discretizationErrorMode=='analytical':
#
# [EOC_l2, EOC_h1, EOC_energy,L2_fine , L2_coarse , H1_fine , H1_coarse, constraintError] \
# = computeError_analytical(numLevels,domainDim,targetDim,order,interpolationMethod,discretizationErrorMode,maxTrustRegionSteps,energyList,ndofList,stepList,initial_energyList,referenceSolution,numReferenceLevels)
#
# else:
# EOC_l2 = "-"
# EOC_h1 = "-"
# EOC_energy = "-"
# L2_fine = "-"
# L2_coarse = "-"
# H1_fine = "-"
# H1_coarse = "-"
# constraintError = "-"
#
# # if maxTrustRegionSteps > 0 :
# # print('energyList', energyList)
# # # currentEnergy = energyList[numLevels-1]
# # currentEnergy = energyList[-1] # better like this?
# # print('currentEnergy', currentEnergy)
# #
# # # currentDofs = ndofList[numLevels-1]
# # currentDofs = ndofList[-1]
# # print('currentDofs', currentDofs)
# #
# # # currentSteps = stepList[numLevels-1]
# # currentSteps = stepList[-1]
# # print('currentSteps', currentSteps)
# #
# # print('initial_energyList', initial_energyList)
# # currentInitialEnergy = initial_energyList[-1]
# # print('currentInitialEnergy', currentInitialEnergy)
# #
# # else:
# # currentEnergy = "-"
# # currentDofs = "-"
# # currentSteps = "-"
# # currentInitialEnergy = "-"
#
# # x.add_row([numLevels, currentDofs , currentSteps, currentInitialEnergy, currentEnergy, "-", "-" , "-" , "-", "-", "-", "-"])
# # print('Energy List:', energyList)
#
#
# if maxTrustRegionSteps > 0 :
# print('energyList', energyList)
# # currentEnergy = energyList[numLevels-1]
# currentEnergy = energyList[-1] # better like this?
# print('currentEnergy', currentEnergy)
#
# # currentDofs = ndofList[numLevels-1]
# currentDofs = ndofList[-1]
# print('currentDofs', currentDofs)
#
# # currentSteps = stepList[numLevels-1]
# currentSteps = stepList[-1]
# print('currentSteps', currentSteps)
#
# print('initial_energyList', initial_energyList)
# currentInitialEnergy = initial_energyList[-1]
# print('currentInitialEnergy', currentInitialEnergy)
#
#
#
#
# else:
# currentEnergy = "-"
# currentDofs = "-"
# currentSteps = "-"
# currentInitialEnergy = "-"
#
# print('timeList:', timeList)
# currentTime = timeList[-1]
# print('currentSteps:', currentSteps)
#
# x.add_row([numLevels, currentDofs, currentSteps, currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy,L2_fine , L2_coarse , H1_fine , H1_coarse , currentTime])
# rows.append([numLevels, currentDofs, str(currentSteps), currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy , currentTime])
#
# ## Add extra column:
# # print(*EOC_l2_list)
# # x.add_column('EOC', EOC_l2_list)
#
# # Write Table to Text-File:
# tablefile = open("EOC-table_R"+str(domainDim)+"intoR"+ str(targetDim)+ "_deg" + str(kappa) + "_" + interpolationMethod+ "_order" +str(order) + "tol" + str(tolerance) + ".txt", "w")
# tablefile.write(str(x))
# tablefile.write('\n')
# tablefile.write(str(tabulate(rows, headers='firstrow', tablefmt='latex')))
#
# tablefile.close()
# # print Table
# print(x)
#
#
#
# #--- Try Tabulate:
#
# print(tabulate(rows, headers='firstrow', tablefmt='latex'))
# ########## end of kappa-loop ##########################
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment