From fc9c4f4aa8711f13e24dbe87f9b000c85884c929 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 31 Aug 2022 15:49:40 +0200
Subject: [PATCH] Add new material setup routine
---
dune/microstructure/CorrectorComputer.hh | 8 +-
dune/microstructure/materialDefinitions.hh | 49 +++++
dune/microstructure/prestrainedMaterial.hh | 229 ++++++++++++++++++-
geometries/material.py | 243 ++++-----------------
inputs/cellsolver.parset | 13 +-
src/Cell-Problem.cc | 4 +-
6 files changed, 328 insertions(+), 218 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index d9edddfd..cad9d1ae 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -944,9 +944,11 @@ public:
const bool verbose = true;
const unsigned int arppp_a_verbosity_level = 2;
const unsigned int pia_verbosity_level = 1;
- MatrixInfo<MatrixCT> matrixInfo(stiffnessMatrix_,verbose,arppp_a_verbosity_level,pia_verbosity_level);
- std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(true) << std::endl;
-
+ if(parameterSet_.get<bool>("print_conditionNumber", false))
+ {
+ MatrixInfo<MatrixCT> matrixInfo(stiffnessMatrix_,verbose,arppp_a_verbosity_level,pia_verbosity_level);
+ std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(true) << std::endl;
+ }
///////////////////////////////////
// --- Choose Solver ---
// 1 : CG-Solver
diff --git a/dune/microstructure/materialDefinitions.hh b/dune/microstructure/materialDefinitions.hh
index 3f13916d..1a03b883 100644
--- a/dune/microstructure/materialDefinitions.hh
+++ b/dune/microstructure/materialDefinitions.hh
@@ -20,6 +20,55 @@ using VectorRT = FieldVector< double, 3>;
+// (30.8.22) DEPRECATED CLASS!
+
+
+
+ ////----------------------------------------------------------------------------------
+// template<class Domain>
+// int indicatorFunction_material_new(const Domain& x)
+// {
+// double theta=0.25;
+// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
+// return 1; //#Phase1
+// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
+// return 2; //#Phase2
+// else
+// return 3; //#Phase3
+// }
+// MatrixRT material_new(const MatrixRT& G, const int& phase)
+// {
+// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
+// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
+// if (phase == 1)
+// return VoigtToMatrix(isotropic(mu[0],lambda[0])* matrixToVoigt(G));
+// // return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1 //Isotrop(mu,lambda)
+// else if (phase == 2)
+// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
+// else
+// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); //#Phase3
+// }
+// MatrixRT prestrain_material_new(const int& phase)
+// {
+// if (phase == 1)
+// return {{1.0, 0.0, 0.0}, //#Phase1
+// {0.0, 1.0, 0.0},
+// {0.0, 0.0, 1.0}
+// };
+// else if (phase == 2)
+// return {{1.0, 0.0, 0.0}, //#Phase2
+// {0.0, 1.0, 0.0},
+// {0.0, 0.0, 1.0}
+// };
+// else
+// return {{0.0, 0.0, 0.0}, //#Phase3
+// {0.0, 0.0, 0.0},
+// {0.0, 0.0, 0.0}
+// };
+// }
+ ////----------------------------------------------------------------------------------
+
+
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 4c9e489d..e158dac7 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -187,9 +187,13 @@ public:
// L3_ = orthotropic(e1,e2,e3,Nspruce_parameters);
- L1_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
- L2_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
- L3_ = transversely_isotropic(e1,e2,e3, {10.7e3,710 ,500,0.51 ,0.31});
+
+ setup("material");
+
+
+ // L1_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
+ // L2_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
+ // L3_ = transversely_isotropic(e1,e2,e3, {10.7e3,710 ,500,0.51 ,0.31});
// L1_ = isotropic(mu[0],lambda[0]);
// L2_ = isotropic(mu[1],lambda[1]);
@@ -244,6 +248,173 @@ public:
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
+
+
+
+ ///////////////////////////////////////////////////////////////////////
+
+// void elasticitySetter //set L1_...
+
+
+FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
+{
+ std::string materialParameters_string;
+ std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
+
+ switch (phase)
+ {
+ case 1:
+ {
+ materialParameters_string = "materialParameters_phase1";
+ break;
+ }
+ case 2:
+ {
+ materialParameters_string = "materialParameters_phase2";
+ break;
+ }
+ case 3:
+ {
+ materialParameters_string= "materialParameters_phase3";
+ break;
+ }
+ default:
+ DUNE_THROW(Exception, "Phase number not implemented.");
+ break;
+ }
+
+ if(phaseType == "isotropic") //TODO SetupPHASE (phase_type)
+ {
+ FieldVector<double,2> materialParameters(0);
+ module.get(materialParameters_string).toC<FieldVector<double,2>>(materialParameters);
+ return isotropic(materialParameters[0],materialParameters[1]);
+ }
+ if(phaseType == "orthotropic") //TODO SetupPHASE (phase_type)
+ {
+ // TODO: Generalize here:
+ //frame vectors
+ VectorRT e1 = {1.0,0.0,0.0};
+ VectorRT e2 = {0.0,1.0,0.0};
+ VectorRT e3 = {0.0,0.0,1.0};
+ FieldVector<double,9> materialParameters(0);
+ module.get(materialParameters_string).toC<FieldVector<double,9>>(materialParameters);
+ return orthotropic(e1,e2,e3,materialParameters);
+ }
+ if(phaseType == "transversely_isotropic") //TODO SetupPHASE (phase_type)
+ {
+ // TODO: Generalize here:
+ //frame vectors
+ VectorRT e1 = {1.0,0.0,0.0};
+ VectorRT e2 = {0.0,1.0,0.0};
+ VectorRT e3 = {0.0,0.0,1.0};
+ FieldVector<double,5> materialParameters(0);
+ module.get(materialParameters_string).toC<FieldVector<double,5>>(materialParameters);
+ return transversely_isotropic(e1,e2,e3,materialParameters);
+ }
+ if(phaseType == "general_anisotropic") //TODO SetupPHASE (phase_type)
+ {
+ // TODO: Generalize here:
+ //frame vectors
+ VectorRT e1 = {1.0,0.0,0.0};
+ VectorRT e2 = {0.0,1.0,0.0};
+ VectorRT e3 = {0.0,0.0,1.0};
+ FieldMatrix<double,6,6> materialParameters(0); //compliance matric
+ module.get(materialParameters_string).toC<FieldMatrix<double,6,6>>(materialParameters);
+ printmatrix(std::cout, materialParameters, "Compliance matrix used:", "--");
+ return general_anisotropic(e1,e2,e3,materialParameters);
+ }
+ else
+ DUNE_THROW(Exception, " Unknown material class for phase ");
+
+}
+
+
+
+
+ //---function that determines elasticity Tensor
+ // void setup(const std::string name, const int Phases) // cant use materialFunctionName_ here!?
+ void setup(const std::string name) // cant use materialFunctionName_ here!?
+ {
+
+ Python::Module module = Python::import(name);
+ indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
+
+ std::cout << "---- setup material... " << std::endl;
+
+ int Phases;
+ module.get("Phases").toC<int>(Phases);
+ std::cout << "Number of Phases: " << Phases << std::endl;
+
+
+ switch (Phases)
+ {
+ case 1:
+ {
+ std::string phase1_type;
+ module.get("phase1_type").toC<std::string>(phase1_type);
+ L1_ = setupPhase(phase1_type, module,1);
+ break;
+ }
+ case 2:
+ {
+ std::string phase1_type;
+ module.get("phase1_type").toC<std::string>(phase1_type);
+ std::string phase2_type;
+ module.get("phase2_type").toC<std::string>(phase2_type);
+
+ L1_ = setupPhase(phase1_type, module,1);
+ L2_ = setupPhase(phase2_type, module,2);
+
+ break;
+ }
+ case 3:
+ {
+ std::string phase1_type;
+ module.get("phase1_type").toC<std::string>(phase1_type);
+ std::string phase2_type;
+ module.get("phase2_type").toC<std::string>(phase2_type);
+ std::string phase3_type;
+ module.get("phase3_type").toC<std::string>(phase3_type);
+ // std::string phase1_type = parameterSet_.get<std::string>("phase1_type", "isotropic");
+ // std::string phase2_type = parameterSet_.get<std::string>("phase2_type", "isotropic");
+ // std::string phase3_type = parameterSet_.get<std::string>("phase3_type", "isotropic");
+
+ L1_ = setupPhase(phase1_type, module,1);
+ L2_ = setupPhase(phase2_type, module,2);
+ L3_ = setupPhase(phase3_type, module,3);
+
+ // if (phase1_type == "isotropic" ) //TODO SetupPHASE (phase_type)
+ // {
+ // FieldVector<double,2> materialParameters_phase1(0);
+ // module.get("materialParameters_phase1").toC<FieldVector<double,2>>(materialParameters_phase1);
+ // L1_ = isotropic(materialParameters_phase1[0],materialParameters_phase1[1]);
+ // }
+ // if (phase2_type == "isotropic" ) //TODO SetupPHASE (phase_type)
+ // {
+ // FieldVector<double,2> materialParameters_phase2(0);
+ // module.get("materialParameters_phase2").toC<FieldVector<double,2>>(materialParameters_phase2);
+ // L2_ = isotropic(materialParameters_phase2[0],materialParameters_phase2[1]);
+ // }
+ // if (phase3_type == "isotropic" ) //TODO SetupPHASE (phase_type)
+ // {
+ // FieldVector<double,2> materialParameters_phase3(0);
+ // module.get("materialParameters_phase3").toC<FieldVector<double,2>>(materialParameters_phase3);
+ // L3_ = isotropic(materialParameters_phase3[0],materialParameters_phase3[1]);
+ // }
+ break;
+ }
+ default:
+ break;
+ }
+ std::cout << "Material setup done." << std::endl;
+ }
+
+ ///////////////////////////////////////////////////////////////////////
+
+
+
+
+
//---function that determines elasticity Tensor
void setupMaterial(const std::string name) // cant use materialFunctionName_ here!?
{
@@ -276,10 +447,9 @@ public:
DUNE_THROW(Exception, "There exists no material with this name ");
}
-//////////// TESTING //////////////////////////////////////
-
+ // --- Helpers
static FieldVector<double,6> MatrixToVoigt(const MatrixRT& G)
{
return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
@@ -291,6 +461,9 @@ public:
}
+ ////////////////////////////////////////////////////////
+ // MATERIAL CLASSES DEFINITIONS:
+ ////////////////////////////////////////////////////////
//--- Definition of isotropic elasticity tensor (in voigt notation)
FieldMatrix<double,6,6> isotropic(const double& mu, const double& lambda)
@@ -303,6 +476,7 @@ public:
{ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 2.0*mu}
};
}
+ //-------------------------------------------------------------------------------
//--- Definition of orthotropic elasticity tensor (in voigt notation)
// - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
@@ -345,6 +519,7 @@ public:
// printmatrix(std::cout, C, "C after .invert()", "--");
return C;
}
+ //-------------------------------------------------------------------------------
//--- Definition of transversely isotropic elasticity tensor (in voigt notation)
@@ -387,10 +562,25 @@ public:
// printmatrix(std::cout, C, "C after .invert()", "--");
return C;
}
+ //-------------------------------------------------------------------------------
-
-
+ // --- Definition of transversely isotropic elasticity tensor (in voigt notation)
+ // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
+ // - Output: compliance Matrix
+ FieldMatrix<double,6,6> general_anisotropic(const VectorRT& framevector1,
+ const VectorRT& framevector2,
+ const VectorRT& framevector3,
+ const FieldMatrix<double,6,6>& C //compliance matrix
+ )
+ {
+ //--- return elasticity tensor
+ auto tmp = C;
+ tmp.invert();
+ // printmatrix(std::cout, C, "C after .invert()", "--");
+ return tmp;
+ }
+ //-------------------------------------------------------------------------------
@@ -434,6 +624,31 @@ public:
// return {{out[0], out[5], out[4]}, {out[5], out[1], out[3]}, {out[4], out[3], out[2]}};
}
+ // MatrixRT prestrain(const int& phase,const Domain& x) const
+ // {
+
+
+ // if (phase == 1)
+ // {
+ // // printmatrix(std::cout, L1_, "L1_", "--");
+ // prestrain1_.mv(G_tmp,out);
+ // }
+ // else if (phase == 2)
+ // {
+ // // printmatrix(std::cout, L2_, "L2_", "--");
+ // L2_.mv(G_tmp,out);
+ // }
+ // else
+ // {
+ // // printmatrix(std::cout, L3_, "L3_", "--");
+ // L3_.mv(G_tmp,out);
+ // }
+
+
+ // }
+
+
+
diff --git a/geometries/material.py b/geometries/material.py
index ff9bcfe2..d1cc5ade 100644
--- a/geometries/material.py
+++ b/geometries/material.py
@@ -5,237 +5,70 @@ import os
import sys
-Phases = 3
-mu_ = [80, 80, 60]
-lambda_ = [80, 80, 25]
-# print('mu_:', mu_)
-# A = [[1, 5, 0], [5,1,0], [5,0,1]]
-#
-# print("sym(A):", sym(A))
-
-
-# dir_path = os.path.dirname(os.path.realpath("/home/klaus/Desktop/Dune-Testing/dune-microstructure/dune/microstructure/matrix_operations.hh"))
-# handle = ctypes.CDLL(dir_path)
-#
-# handle.create2Darray.argtypes = [ctypes.c_int, ctypes.c_double, ctypes.c_double]
-#
-# def create2Darray(nside, mx, my):
-# return handle.create2Darray(nside, mx, my)
-
-
-#Indicator function that determines both phases
-# x[0] : y1-component
-# x[1] : y2-component
-# x[2] : x3-component
-# --- replace with your definition of indicatorFunction:
-###############
-# Wood
-###############
-# def f(x):
-# theta=0.25
-# # mu_ = [3, 5, 6]
-# # lambda_ = [9, 7, 8]
-# # mu_ = 3 5 6
-# # lambda_ = 9 7 8
-
-# if ((abs(x[0]) < theta/2) and x[2]<0.25):
-# return [mu_[0], lambda_[0]] #latewood
-# # return 5 #latewood
-# elif ((abs(x[0]) > theta/2) and x[2]<0.25):
-# return [mu_[1], lambda_[1]] #latewood
-# # return 2
-# else :
-# return [mu_[2],lambda_[2]] #latewood #Phase3
-# # return 1
-
def indicatorFunction(x):
theta=0.25
factor=1
- if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+ if (x[0] <-0.5+theta and x[2]<-0.5+theta):
return 1 #Phase1
- elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+ elif (x[1]<-0.5+theta and x[2]>0.5-theta):
return 2 #Phase2
else :
return 0 #Phase3
-def L1(G):
- return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
-
-def L2(G):
- return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
-
-def L3(G):
- return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
-
-
-# TEST
-
-# def L1(G):
-# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) )) #Phase1
-
-# def L2(G):
-# return Add(smult(2.0 * mu_[1], sym(G)),smult(lambda_[1] ,smult(trace(sym(G)),Id()) )) #Phase1
-
-# def L3(G):
-# return Add(smult(2.0 * mu_[2], sym(G)),smult(lambda_[2] ,smult(trace(sym(G)),Id()) )) #Phase1
-
-
-#Workaround
-def muValue(x):
- return mu_
+########### Options for material phases: #################################
+# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic"
+#########################################################################
+## Notation (materialParameters_phase):
+# isotropic (Lame parameters): [mu , lambda]
+# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic : [E1,E2,G12,nu12,nu23]
+#
+# general_anisotropic : full compliance matrix C
+######################################################################
-def lambdaValue(x):
- return lambda_
+# --- Number of material phases
+Phases=3
+##--- Define different material phases:
+phase1_type="orthotropic"
+materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37] # walnut_parameters (values for compliance matrix)
+# materialParameters_phase1 = [10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31] # Nspruce_parameters (values for compliance matrix)
+phase2_type="transversely_isotropic"
+materialParameters_phase2 = [11.2e3,1190,960,0.63 ,0.37]
+# phase1_type="isotropic"
+# materialParameters_phase1 = [80, 80]
+#
+# phase2_type="isotropic"
+# materialParameters_phase2 = [80, 80]
-###############
-# Cross
-###############
-# def f(x):
-# theta=0.25
-# factor=1
-# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
-# return 1 #Phase1
-# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
-# return 2 #Phase2
-# else :
-# return 0 #Phase3
+# phase3_type="isotropic"
+# materialParameters_phase3 = [60, 25]
+# for general anisotopic material the compliance matrix is required:
+phase3_type="general_anisotropic"
+materialParameters_phase3 = np.array([[1.0, 0.0, 0.0, 0.0 , 0.0, 0.0],
+ [0.0, 1.0, 0.0, 0.0 , 0.0, 0.0],
+ [0.0, 0.0, 1.0, 0.0 , 0.0, 0.0],
+ [0.0, 0.0, 0.0, math.sqrt(2), 0.0, 0.0],
+ [0.0, 0.0, 0.0, 0.0 , 1.0, 0.0],
+ [0.0, 0.0, 0.0, 0.0 , 0.0, 1.0]])
-# def f(x):
-# # --- replace with your definition of indicatorFunction:
-# if (abs(x[0]) > 0.25):
-# return 1 #Phase1
-# else :
-# return 0 #Phase2
-def b1(x):
+def prestrain_phase1(x):
return [[1, 0, 0], [0,1,0], [0,0,1]]
-def b2(x):
+def prestrain_phase2(x):
return [[1, 0, 0], [0,1,0], [0,0,1]]
-def b3(x):
+def prestrain_phase3(x):
return [[0, 0, 0], [0,0,0], [0,0,0]]
-
-# mu=80 80 60
-
-# lambda=80 80 25
-
-
-# --- elasticity tensor
-# def L(G,x):
-# def L(G):
-# # input: symmetric matrix G, position x
-# # output: symmetric matrix LG
-# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
-
-
-
-
-
-
-
-
-
-# --- elasticity tensor
-def L(G,x):
- # input: symmetric matrix G, position x
- # output: symmetric matrix LG
- theta=0.25
- if (x[0] <-1/2+theta and x[2]<-1/2+theta):
- return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
- elif (x[1]< -1/2+theta and x[2]>1/2-theta):
- return 2.0 * mu_[1] * sym(G) + lambda_[1] * trace(sym(G)) * Id() #Phase2
- else :
- return 2.0 * mu_[2] * sym(G) + lambda_[2] * trace(sym(G)) * Id() #Phase3
-# # # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
-
-
-# def L(G,x):
-# # input: symmetric matrix G, position x
-# # output: symmetric matrix LG
-# theta=0.25
-# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
-# # return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
-# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) ))
-# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
-# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) )) #Phase2
-# else :
-# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) )) #Phase3
-# # return [[0, 0, 0], [0,0,0], [0,0,0]] #Phase3
-
-
-##TEST
-# def L(G,x):
-# # input: symmetric matrix G, position x
-# # output: symmetric matrix LG
-# theta=0.25
-# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
-# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
-# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
-# return [[x[0], 1, x[0]], [1, 1, 1],[x[0], x[0], x[0]]]
-# else :
-# return [[0, x[2], x[2]], [0,x[2],0], [0,0,0]]
-
-
-##TEST
-# def L(G,x):
-# # input: symmetric matrix G, position x
-# # output: symmetric matrix LG
-# theta=0.25
-# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
-# return sym([[1, 1, 1], [1, 1, 1],[1, 1, 1]])
-# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
-# return sym([[x[0], 1, x[0]], [1, 1, 1],[x[0], x[0], x[0]]])
-# else :
-# return sym([[0, x[2], x[2]], [0,x[2],0], [0,0,0]])
-
-
-
-# # small speedup..
-# def L(G,x):
-# # input: symmetric matrix G, position x
-# # output: symmetric matrix LG
-# theta=0.25
-# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
-# return mu_[0] * (np.array(G).transpose() + np.array(G)) + lambda_[0] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3) #Phase1
-# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
-# return mu_[1] * (np.array(G).transpose() + np.array(G)) + lambda_[1] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3) #Phase2
-# else :
-# return mu_[2] * (np.array(G).transpose() + np.array(G)) + lambda_[2] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3) #Phase3
-# # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
-
-
-
-
-
-
-# def H(G,x):
-# # input: symmetric matrix G, position x
-# # output: symmetric matrix LG
-# if (abs(x[0]) > 0.25):
-# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
-# else:
-# return [[0, 0, 0], [0,0,0], [0,0,0]]
-
-
-def H(G,x):
- # input: symmetric matrix G, position x
- # output: symmetric matrix LG
- if (abs(x[0]) > 0.25):
- return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
- else:
- return [[0, 0, 0], [0,0,0], [0,0,0]]
-
-# 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 88cfed92..c88beabb 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -16,6 +16,8 @@ outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
print_debug = true #(default=false)
+
+
#############################################
# Grid parameters
#############################################
@@ -44,6 +46,11 @@ gamma=1.0
+#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
+geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
+
+
+
#############################################
# Assembly options
@@ -59,7 +66,7 @@ gamma=1.0
#############################################
# Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
#############################################
-Solvertype = 3
+Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
Solver_verbosity = 0 #(default = 2) degree of information for solver output
@@ -86,5 +93,9 @@ write_checkOrthogonality = true
#write_corrector_phi2 = true
#write_corrector_phi3 = true
+
+# --- Print Condition number of matrix (can be expensive):
+#print_conditionNumber= true #(default=false)
+
# --- write effective quantities to Matlab-folder for symbolic minimization:
write_toMATLAB = true # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index ea11074b..7d7f2c3b 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -143,14 +143,14 @@ int main(int argc, char *argv[])
// parameterSet.report(log); // short Alternativ
//--- Get Path for Material/Geometry functions
- // std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath",);
+ std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
//--- Start Python interpreter
Python::start();
Python::Reference main = Python::import("__main__");
Python::run("import math");
Python::runStream()
<< std::endl << "import sys"
- // << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
+ << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
<< std::endl;
--
GitLab