diff --git a/inputs/computeMuGamma.parset b/inputs/computeMuGamma.parset
new file mode 100644
index 0000000000000000000000000000000000000000..6971b9673e304e388ddc8bdce91ba5790d046064
--- /dev/null
+++ b/inputs/computeMuGamma.parset
@@ -0,0 +1,102 @@
+
+#path for logfile
+#outputPath = "../../outputs/output.txt"
+#outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt"
+#outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs"
+
+#############################################
+# Debug Output
+#############################################
+
+#print_debug = true #default = false
+
+
+
+
+
+#############################################
+# Grid parameters
+#############################################
+
+
+nElements = 32 32
+
+gamma=0.5
+#############################################
+# Material parameters
+#############################################
+
+write_materialFunctions = true # VTK mu-functions , lambda-functions
+
+beta = 2.0 # ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case
+
+mu1=10.0
+#mu1=1000.0
+
+#lambda1=10.0
+#lambda1 = 20.0
+#lambda1 = 20.0
+#lambda1 = 5.0
+
+#### material_implementation("analytical_Example") ?
+
+material_prestrain_imp= "analytical_Example"
+#material_prestrain_imp ="isotropic_bilayer"
+
+
+#material_prestrain_imp= "circle_fiber" #TEST
+#material_prestrain_imp= "square_fiber" #TEST
+
+#############################################
+# Prestrain parameters
+#############################################
+
+write_prestrainFunctions = true # VTK norm of B ,
+
+rho1 = 1.0
+alpha = 2.0 # ratio between prestrain parameters rho1 & rho2
+
+#theta = 0.3 # volume fraction #default = 1.0/4.0
+#theta = 0.25 # volume fraction
+#theta = 0.75 # volume fraction
+
+
+
+------------Matrix Material --------------
+#material_prestrain_imp ="matrix_material_circles"
+#material_prestrain_imp ="matrix_material_squares"
+
+nF = 8 #number of Fibers (in each Layer)
+#rF = 0.05 #Fiber radius max-fiber-radius = (width/(2.0*nF)
+------------------------------------------
+
+width = 1.0
+height = 1.0
+
+
+
+
+
+#############################################
+# Solver Type
+#############################################
+# 1: CG-Solver # 2: GMRES # 3: QR
+
+Solvertype = 1
+
+
+
+
+#############################################
+# Define Analytic Solutions
+#############################################
+
+#b1 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2)
+
+
+
+
+
+
+
+
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index bf2c1f8b7b83b1f1053cbfd966cbc925b08ee502..2657bbc0221955e6e04bb8b7982b9dd597f54d5e 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -885,6 +885,7 @@ int main(int argc, char *argv[])
// Get Material Parameters
///////////////////////////////////
std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
+ log << "material_prestrain used: "<< imp << std::endl;
double beta = parameterSet.get<double>("beta",2.0);
double mu1 = parameterSet.get<double>("mu1",10.0);;
double mu2 = beta*mu1;
@@ -1475,8 +1476,12 @@ int main(int argc, char *argv[])
double b2_eff_ana = (-3.0/2.0)*(theta*mu2)/(mu1*(1-theta)+mu2*theta);
double b3_eff_ana = 0.0;
- double q1_ana = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
- double q2_ana = ((1.0-theta)*mu1+theta*mu2)/6.0;
+// double q1_ana = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+// double q2_ana = ((1.0-theta)*mu1+theta*mu2)/6.0;
+ double q1_ana = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0); // 1/6 * harmonicMean
+ double q2_ana = mu1*((1-theta)+theta*beta) *(1.0/6.0); // 1/6 * arithmeticMean
+
+
@@ -1603,7 +1608,7 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
- std::string VTKOutputName = "CellProblem-result";
+ std::string vtkOutputName = outputPath + "/CellProblem-result";
VTKWriter<GridView> vtkWriter(gridView_CE);
@@ -1616,10 +1621,10 @@ int main(int argc, char *argv[])
vtkWriter.addVertexData(
correctorFunction_3,
VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
- vtkWriter.write( VTKOutputName + "-level"+ std::to_string(level));
-// vtkWriter.pwrite( VTKOutputName + "-level"+ std::to_string(level), outputPath, ""); // TEST Write to folder "/outputs"
-// vtkWriter.pwrite( VTKOutputName + "-level"+ std::to_string(level), outputPath, "", VTK::OutputType::ascii, 0, 0 );
- std::cout << "wrote data to file: " + VTKOutputName + "-level" + std::to_string(level) << std::endl;
+ vtkWriter.write( vtkOutputName + "-level"+ std::to_string(level));
+// vtkWriter.pwrite( vtkOutputName + "-level"+ std::to_string(level), outputPath, ""); // TEST Write to folder "/outputs"
+// vtkWriter.pwrite( vtkOutputName + "-level"+ std::to_string(level), outputPath, "", VTK::OutputType::ascii, 0, 0 );
+ std::cout << "wrote data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
@@ -1665,8 +1670,8 @@ int main(int argc, char *argv[])
lambda_DGBF_P0,
VTK::FieldInfo("lambda_P0", VTK::FieldInfo::Type::scalar, 1));
- MaterialVtkWriter.write("MaterialFunctions-level"+ std::to_string(level) );
- std::cout << "wrote data to file: MaterialFunctions-level" + std::to_string(level) << std::endl;
+ MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
+ std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;
}
diff --git a/src/Compute_MuGamma.cc b/src/Compute_MuGamma.cc
index 3cef8ce93319663790cb25915e468b82f4450186..17f553785d0eaa573a56212f1bc7dd1c23d5a314 100644
--- a/src/Compute_MuGamma.cc
+++ b/src/Compute_MuGamma.cc
@@ -891,24 +891,22 @@ int main(int argc, char *argv[])
double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
-// double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);
-// double am = mu1*((1-theta)+theta*beta) *(1.0/6.0);
+ double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);
+ double am = mu1*((1-theta)+theta*beta) *(1.0/6.0);
std::cout << "q1 : " << q1 << std::endl;
std::cout << "q2 : " << q2 << std::endl;
std::cout << "q3 should be between q1 and q2" << std::endl;
-// std::cout << "hm : " << hm << std::endl;
-// std::cout << "am : " << am << std::endl;
+ std::cout << "hm : " << hm << std::endl;
+ std::cout << "am : " << am << std::endl;
if(mu_gamma > q2)
{
- std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
+ std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
}
-
- std::cout << "beta : " << beta << std::endl;
std::cout << "beta : " << beta << std::endl;
std::cout << "theta : " << theta << std::endl;
std::cout << "Gamma : " << gamma << std::endl;
@@ -918,6 +916,8 @@ int main(int argc, char *argv[])
// Output result
+ std::string VTKOutputName = outputPath + "/Compute_MuGamma-Result";
+
VTKWriter<GridView> vtkWriter(gridView);
// vtkWriter.addVertexData(x, "solution"); //--- Anpassen für P2
vtkWriter.addVertexData(
@@ -925,7 +925,8 @@ int main(int argc, char *argv[])
VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
// VTK::FieldInfo("solution", VTK::FieldInfo::Type::vector, dim));
- vtkWriter.write("Dune-Poisson-Result");
+ vtkWriter.write( VTKOutputName );
+ std::cout << "wrote data to file: " + VTKOutputName << std::endl;
@@ -946,7 +947,7 @@ int main(int argc, char *argv[])
mu_DGBF_P0,
VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));
- MaterialVtkWriter.write("MaterialFunctions-level" );
- std::cout << "wrote data to file: MaterialFunctions" << std::endl;
+ MaterialVtkWriter.write(outputPath + "/MaterialFunctions" );
+ std::cout << "wrote data to file:" + outputPath + "/MaterialFunctions" << std::endl;
}
diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py
index 796aee1e6aef428819a67da97e1f0eff853bbbec..5fec95a9955f772746d14c954c199988f3e384b1 100644
--- a/src/PhaseDiagram.py
+++ b/src/PhaseDiagram.py
@@ -9,83 +9,185 @@ import re
import matlab.engine
import sys
from ClassifyMin import *
+# from CellScript import *
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
from vtk.util import numpy_support
from pyevtk.hl import gridToVTK
+import time
# print(sys.executable)
-# from subprocess import Popen, PIPE
-# Not Needed:
-# def ndm(*args):
-# return [x[(None,)*i+(slice(None),)+(None,)*(len(args)-i-1)] for i, x in enumerate(args)]
+# --------------------------------------------------------------------
+# START :
+# INPUT (Parameters): alpha, beta, theta, gamma, mu1, rho1
+#
+# -Option 1 : (Case lambda = 0 => q12 = 0)
+# compute q1,q2,b1,b2 from Formula
+# Option 1.1 :
+# set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
+# Option 1.2 :
+# compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
+# -Option 2 :
+# compute Q_hom & B_eff by running 'Cell-Problem'
+#
+# -> CLASSIFY ...
+#
+# OUTPUT: Minimizer G, angle , type, curvature
+# -----------------------------------------------------------------------
+
+
+# unabhängig von alpha...
+def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+ # ------------------------------------ get mu_gamma ------------------------------
+ # ---Scenario 1.1: extreme regimes
+ if gamma == '0':
+ print('extreme regime: gamma = 0')
+ mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+ print("mu_gamma:", mu_gamma)
+ elif gamma == 'infinity':
+ print('extreme regime: gamma = infinity')
+ mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta) # = q1
+ print("mu_gamma:", mu_gamma)
+ else:
+ # --- Scenario 1.2: compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+ print("Run computeMuGamma for Gamma = ", gamma)
+ # print('gamma='+str(gamma))
+ with open(InputFilePath, 'r') as file:
+ filedata = file.read()
+ filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+ # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+ filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+ filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+ filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+ filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+ f = open(InputFilePath,'w')
+ f.write(filedata)
+ f.close()
+ # --- Run Cell-Problem
+ t = time.time()
+ # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+ # capture_output=True, text=True)
+ # --- Run Cell-Problem_muGama -> faster
+ # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+ # capture_output=True, text=True)
+ # --- Run Compute_muGamma (2D Problem much much faster)
+ subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+ capture_output=True, text=True)
+ print('elapsed time:', time.time() - t)
+
+ #Extract mu_gamma from Output-File TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+ with open(OutputFilePath, 'r') as file:
+ output = file.read()
+ tmp = re.search(r'(?m)^mu_gamma=.*',output).group() # Not necessary for Intention of Program t output Minimizer etc.....
+ s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+ mu_gamma = float(s[0])
+ # print("mu_gamma:", mu_gammaValue)
+ # --------------------------------------------------------------------------------------
+ return mu_gamma
+
+
+
+
+
+# --- SETUPS PATHS
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
-# ---------------------------------------------- Main ---------------------s
-# --- Input Parameters ----
-mu_1 = 1.0
-rho_1 = 1.0
-alpha = 9.0
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+
+# -------------------------- Input Parameters --------------------
+mu1 = 10.0 # TODO : here must be the same values as in the Parset for computeMuGamma
+rho1 = 1.0
+alpha = 2.0
beta = 2.0
-theta = 1.0/8.0
-# print('Type of theta', type(theta))
+theta = 1.0/4.0
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+gamma = '0'
+gamma = 'infinity'
+gamma = 0.5
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
+# muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
+
+print('Test MuGamma:', muGamma)
# ------- Options --------
# print_Cases = True
# print_Output = True
+
# make_3D_plot = True
-# make_2D_plot = False
# make_3D_PhaseDiagram = True
+# make_2D_plot = False
# make_2D_PhaseDiagram = False
-
make_3D_plot = False
-make_2D_plot = True
make_3D_PhaseDiagram = False
+make_2D_plot = True
make_2D_PhaseDiagram = True
-# --- Define q1, q2 , mu_gamma, q12 ---
-# 1. read from Cell-Output
-# 2. define values from analytic formulas (expect for mu_gamma)
-q1 = harmonicMean(mu_1, beta, theta)
-q2 = arithmeticMean(mu_1, beta, theta)
+# --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
+# q1 = harmonicMean(mu1, beta, theta)
+# q2 = arithmeticMean(mu1, beta, theta)
+# --- Set q12
+# q12 = 0.0 # (analytical example) # TEST / TODO read from Cell-Output
-# TEST / TODO read from Cell-Output
-q12 = 0.0 # (analytical example)
-# TODO read from Cell-Output
-# --- Set mu_gamma to value
-mu_gamma = q1 # TODO read from Cell-Output
-b1 = prestrain_b1(rho_1, beta, alpha, theta)
-b2 = prestrain_b2(rho_1, beta, alpha, theta)
-print('---- Input parameters: -----')
-print('mu_1: ', mu_1)
-print('rho_1: ', rho_1)
-print('alpha: ', alpha)
-print('beta: ', beta)
-print('theta: ', theta)
-print("q1: ", q1)
-print("q2: ", q2)
-print("mu_gamma: ", mu_gamma)
-print("q12: ", q12)
-print("b1: ", b1)
-print("b2: ", b2)
-print('----------------------------')
+# b1 = prestrain_b1(rho1, beta, alpha, theta)
+# b2 = prestrain_b2(rho1, beta, alpha, theta)
+#
+# print('---- Input parameters: -----')
+# print('mu1: ', mu1)
+# print('rho1: ', rho1)
+# print('alpha: ', alpha)
+# print('beta: ', beta)
+# print('theta: ', theta)
+# print("q1: ", q1)
+# print("q2: ", q2)
+# print("mu_gamma: ", mu_gamma)
+# print("q12: ", q12)
+# print("b1: ", b1)
+# print("b2: ", b2)
+# print('----------------------------')
# print("machine epsilon", sys.float_info.epsilon)
-
# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
-
# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
# print("Test", Test)
# ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
-SamplePoints_3D = 100 # Number of sample points in each direction
-SamplePoints_2D = 800 # Number of sample points in each direction
+SamplePoints_3D = 10 # Number of sample points in each direction
+SamplePoints_2D = 10 # Number of sample points in each direction
+SamplePoints_3D = 40 # Number of sample points in each direction
+SamplePoints_2D = 3 # Number of sample points in each direction
+
+
if make_3D_PhaseDiagram:
alphas_ = np.linspace(-20, 20, SamplePoints_3D)
@@ -94,14 +196,19 @@ if make_3D_PhaseDiagram:
# print('type of alphas', type(alphas_))
# print('Test:', type(np.array([mu_gamma])) )
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
- classifyMin_geoVec = np.vectorize(classifyMin_geo)
- G, angles, Types, curvature = classifyMin_geoVec(alphas,betas,thetas,mu_gamma,q12, mu_1, rho_1)
+ classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+ # Get MuGamma values ...
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1)
+ # Classify Minimizers....
+ G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas,mu_gamma, mu1, rho1) # Sets q12 to zero!!!
# print('size of G:', G.shape)
# print('G:', G)
- # Out = classifyMin_geoVec(alphas,betas,thetas)
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
# T = Out[2]
# --- Write to VTK
- VTKOutputName = "./PhaseDiagram3D"
+ VTKOutputName = "outputs/PhaseDiagram3D"
gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
print('Written to VTK-File:', VTKOutputName )
@@ -114,14 +221,20 @@ if make_2D_PhaseDiagram:
# print('type of alphas', type(alphas_))
# print('Test:', type(np.array([mu_gamma])) )
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
- classifyMin_geoVec = np.vectorize(classifyMin_geo)
- G, angles, Types, curvature = classifyMin_geoVec(alphas,betas,thetas,mu_gamma,q12, mu_1, rho_1)
+
+
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+ muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1)
+ G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
# print('size of G:', G.shape)
# print('G:', G)
- # Out = classifyMin_geoVec(alphas,betas,thetas)
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
# T = Out[2]
# --- Write to VTK
- VTKOutputName = "./PhaseDiagram2D"
+ # VTKOutputName = + path + "./PhaseDiagram2DNEW"
+ VTKOutputName = "outputs/PhaseDiagram2D"
gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
print('Written to VTK-File:', VTKOutputName )