diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh index 0d56aade9b06071c96bb9669a87a2de617fe456d..9329f1634fda83af82f0b28c2c896b99aa859505 100644 --- a/dune/microstructure/EffectiveQuantitiesComputer.hh +++ b/dune/microstructure/EffectiveQuantitiesComputer.hh @@ -60,7 +60,7 @@ public: { // prestrain_ = material_.getPrestrainFunction(); - std::cout << "starting effective quantities computation..." << std::endl; + } @@ -79,6 +79,8 @@ public: void computeEffectiveQuantities() { + std::cout << "starting effective quantities computation..." << std::endl; + // Get everything.. better TODO: with Inheritance? // auto test = correctorComputer_.getLoad_alpha1(); // auto phiContainer = correctorComputer_.getPhicontainer(); @@ -179,6 +181,7 @@ public: elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ??? if (b==0) { + // std::cout << "WORKS TILL HERE" << std::endl; // auto prestrainPhaseValue_old = prestrain_(localIndicatorFunction(quadPos)); auto prestrainPhaseValue = material_.prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos)); // if (localIndicatorFunction(quadPos) != 3) @@ -291,7 +294,8 @@ public: const auto& quadPos = quadPoint.position(); const double integrationElement = geometry.integrationElement(quadPos); - double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos))); + double energyDensity= scalarProduct(material_.applyElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos))); + // double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos))); // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos)); elementEnergy += energyDensity * quadPoint.weight() * integrationElement; } diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh index 4a020ed988ec52eedfe97a90bfe2e37a6464a298..8e515fe1822b8cc2ea07a90ef36a818e89eeca05 100644 --- a/dune/microstructure/prestrainedMaterial.hh +++ b/dune/microstructure/prestrainedMaterial.hh @@ -30,9 +30,6 @@ using std::make_shared; - - - // // ---------------------------------------------------------------------------------- // template<class Domain> // int indicatorFunction_material(const Domain& x) @@ -635,8 +632,11 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m return prestrain1_(x); else if (phase == 2) return prestrain2_(x); - else + else if (phase ==3) return prestrain3_(x); + else + DUNE_THROW(Exception, "Phase number not implemented."); + } diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset index 9ce849f0fcccab1c7792360a57e7895b916cecf2..902e0b882f4972fc7714e0c7f97c122438b7ffc5 100644 --- a/inputs/cellsolver.parset +++ b/inputs/cellsolver.parset @@ -26,7 +26,7 @@ outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs ## {start,finish} computes on all grid from 2^(start) to 2^finish refinement #---------------------------------------------------- -numLevels= 2 3 # computes all levels from first to second entry +numLevels= 2 2 # computes all levels from first to second entry ############################################# @@ -34,12 +34,14 @@ numLevels= 2 3 # computes all levels from first to second entry ############################################# # --- Choose material definition: +#materialFunction = "material_test" +#materialFunction = "material_neukamm" #materialFunction = "two_phase_material_1" #materialFunction = "two_phase_material_2" -#materialFunction = "material_neukamm" -materialFunction = "material" -#materialFunction = "material_2" -#materialFunction = "homogeneous" +#materialFunction = "two_phase_material_3" +#materialFunction = "material_orthotropic" +materialFunction = "homogeneous" + # --- Choose scale ratio gamma: gamma=1.0 @@ -47,8 +49,8 @@ gamma=1.0 -#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries -geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries +#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/materials +geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/materials diff --git a/materials/__pycache__/material_neukamm.cpython-38.pyc b/materials/__pycache__/material_neukamm.cpython-38.pyc index f077d3db1ae6d8c15bcd7f2a1addeacdd7531020..9d7421aeea4ce841358b2a6f70980f04d0e1084f 100644 Binary files a/materials/__pycache__/material_neukamm.cpython-38.pyc and b/materials/__pycache__/material_neukamm.cpython-38.pyc differ diff --git a/materials/__pycache__/material_test.cpython-38.pyc b/materials/__pycache__/material_test.cpython-38.pyc index 1dfcaeea6e9e18505d515bf97694734c7d04784c..5185825128d19bce06d794b7260599da858ed7ed 100644 Binary files a/materials/__pycache__/material_test.cpython-38.pyc and b/materials/__pycache__/material_test.cpython-38.pyc differ diff --git a/materials/__pycache__/two_phase_material_1.cpython-38.pyc b/materials/__pycache__/two_phase_material_1.cpython-38.pyc index c9763e8459ecb9c81f61150564e7b902e4558747..78b1dd2bd3091a8f4709ecd09742ae4c3d6ac0a7 100644 Binary files a/materials/__pycache__/two_phase_material_1.cpython-38.pyc and b/materials/__pycache__/two_phase_material_1.cpython-38.pyc differ diff --git a/materials/__pycache__/two_phase_material_2.cpython-38.pyc b/materials/__pycache__/two_phase_material_2.cpython-38.pyc index de9026c11a1592fe0b109938ac6617277d4588af..86b1c59d67e8cc248fd7ae49f0d16bdafe0fc4fc 100644 Binary files a/materials/__pycache__/two_phase_material_2.cpython-38.pyc and b/materials/__pycache__/two_phase_material_2.cpython-38.pyc differ diff --git a/materials/homogeneous.py b/materials/homogeneous.py new file mode 100644 index 0000000000000000000000000000000000000000..7369270db99d37ad7d6d563d6afcb77a0dff4c58 --- /dev/null +++ b/materials/homogeneous.py @@ -0,0 +1,35 @@ +import math + + +#Indicator function that determines both phases +# x[0] : x-component +# x[1] : y-component +# x[2] : z-component +def indicatorFunction(x): + # --- replace with your definition of indicatorFunction: + return 1 #Phase1 + + + +########### Options for material phases: ################################# +# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic" +######################################################################### +## Notation - Parameter input : +# 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 +###################################################################### +# --- Number of material phases +Phases=1 +#--- Define different material phases: +#- PHASE 1 +phase1_type="isotropic" +materialParameters_phase1 = [80, 80] + + + +#--- define prestrain function for each phase +# (also works with non-constant values) +def prestrain_phase1(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] diff --git a/materials/material_neukamm.py b/materials/material_neukamm.py index 42440c6d328bdebf505df8075bcf38b8e6fa512c..85c128f723919481aa44d4ada2802ae44a46aeae 100644 --- a/materials/material_neukamm.py +++ b/materials/material_neukamm.py @@ -1,11 +1,30 @@ import math +from python_matrix_operations import * +import ctypes +import os +import sys +#--------------------------------------------------------------- -#Indicator function that determines both phases + +#--- define indicator function for material phases # x[0] : y1-component # x[1] : y2-component # x[2] : x3-component -# --- replace with your definition of indicatorFunction: +#To indicate phases return either : 1 / 2 / 3 +############### +# Cross +############### +def indicatorFunction(x): + theta=0.25 + factor=1 + if (x[0] <-0.5+theta and x[2]<-0.5+theta): + return 1 #Phase1 + elif (x[1]<-0.5+theta and x[2]>0.5-theta): + return 2 #Phase2 + else : + return 3 #Phase3 + ############### # Wood ############### @@ -18,32 +37,40 @@ import math # else : # return 0 #Phase3 -# def b1(x): -# return [[.5, 0, 0], [0,1,0], [0,0,0]] -# def b2(x): -# return [[.4, 0, 0], [0,.4,0], [0,0,0]] +########### Options for material phases: ################################# +# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic" +######################################################################### +## Notation - Parameter input : +# 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 +###################################################################### +# --- Number of material phases +Phases=3 +#--- Define different material phases: +#- PHASE 1 +phase1_type="isotropic" +materialParameters_phase1 = [80, 80] + +#- PHASE 2 +phase2_type="isotropic" +materialParameters_phase2 = [80, 80] + +#- PHASE 3 +phase3_type="isotropic" +materialParameters_phase3 = [60, 25] -# def b3(x): -# return [[0, 0, 0], [0,0,0], [0,0,0]] -############### -# 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 -def b1(x): +#--- define prestrain function for each phase +# (also works with non-constant values) +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]] + # return [[x[0],0 ,0 ], [0,0,x[1]], [0,0,x[2]]] diff --git a/materials/material_neukamm_old.py b/materials/material_neukamm_old.py new file mode 100644 index 0000000000000000000000000000000000000000..790e3a436fa63b2071fd95815a0a0bf22ab84c6a --- /dev/null +++ b/materials/material_neukamm_old.py @@ -0,0 +1,50 @@ +import math + +# DEPRECATED!!!! just for reference + +#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 + # if ((abs(x[0]) < theta/2) and x[2]<0.25): + # return 1 #latewood + # elif ((abs(x[0]) > theta/2) and x[2]<0.25): + # return 2 #earlywood + # else : + # return 0 #Phase3 + +# def b1(x): +# return [[.5, 0, 0], [0,1,0], [0,0,0]] + +# def b2(x): +# return [[.4, 0, 0], [0,.4,0], [0,0,0]] + +# def b3(x): +# return [[0, 0, 0], [0,0,0], [0,0,0]] +############### +# 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 + +def b1(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def b2(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def b3(x): + return [[0, 0, 0], [0,0,0], [0,0,0]] diff --git a/materials/material_orthotropic.py b/materials/material_orthotropic.py new file mode 100644 index 0000000000000000000000000000000000000000..6e68cb6451034044f23c66b8c132a96d138fdf4d --- /dev/null +++ b/materials/material_orthotropic.py @@ -0,0 +1,71 @@ +import math +from python_matrix_operations import * +import ctypes +import os +import sys +#--------------------------------------------------------------- + + + +#--- define indicator function +def indicatorFunction(x): + theta=0.25 + factor=1 + if (x[0] <-0.5+theta and x[2]<-0.5+theta): + return 1 #Phase1 + elif (x[1]<-0.5+theta and x[2]>0.5-theta): + return 2 #Phase2 + else : + return 3 #Phase3 + + +########### Options for material phases: ################################# +# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic" +######################################################################### +## Notation - Parameter input : +# 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 +###################################################################### + +# --- Number of material phases +Phases=3 + + +#--- Define different material phases: + +#- PHASE 1 +phase1_type="orthotropic" +materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37] # walnut parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109] + +#- PHASE 2 +phase2_type="orthotropic" +materialParameters_phase2 = [10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31] # Norway spruce parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109] + +#- PHASE 3 +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]]) + + +#--- define prestrain function for each phase +# (also works with non-constant values) +def prestrain_phase1(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def prestrain_phase2(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def prestrain_phase3(x): + return [[0, 0, 0], [0,0,0], [0,0,0]] + # return [[x[0],0 ,0 ], [0,0,x[1]], [0,0,x[2]]] diff --git a/materials/material_test.py b/materials/material_test.py index fed02ba54e0513084d34d13913c4992eb3fd4ede..7d59bae264b5b49604cb6165fe30fd566fbeea07 100644 --- a/materials/material_test.py +++ b/materials/material_test.py @@ -5,7 +5,7 @@ import os import sys #--------------------------------------------------------------- - +#To indicate phases return either : 1 / 2 / 3 #--- define indicator function def indicatorFunction(x): diff --git a/materials/two_phase_material_1.py b/materials/two_phase_material_1.py index 8c3418f48755232a24c143d57ffb1d3f842746cf..b09bf3c98204d17019b84ad65cfb42dbee6434f6 100644 --- a/materials/two_phase_material_1.py +++ b/materials/two_phase_material_1.py @@ -5,9 +5,40 @@ import math # x[0] : x-component # x[1] : y-component # x[2] : z-component -def f(x): +def indicatorFunction(x): # --- replace with your definition of indicatorFunction: if (abs(x[0]) > 0.25): return 1 #Phase1 else : - return 0 #Phase2 + return 2 #Phase2 + + +########### Options for material phases: ################################# +# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic" +######################################################################### +## Notation - Parameter input : +# 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 +###################################################################### +# --- Number of material phases +Phases=2 +#--- Define different material phases: +#- PHASE 1 +phase1_type="isotropic" +materialParameters_phase1 = [80, 80] + +#- PHASE 2 +phase2_type="isotropic" +materialParameters_phase2 = [80, 80] + + + +#--- define prestrain function for each phase +# (also works with non-constant values) +def prestrain_phase1(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def prestrain_phase2(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] diff --git a/materials/two_phase_material_2.py b/materials/two_phase_material_2.py index b643b490f72d73dd77f9a72bd2d55a1e5fca9101..b26dd94b7b68934ce5b75f6c19607e3607802a0b 100644 --- a/materials/two_phase_material_2.py +++ b/materials/two_phase_material_2.py @@ -5,9 +5,40 @@ import math # x[0] : x-component # x[1] : y-component # x[2] : z-component -def f(x): +def indicatorFunction(x): # --- replace with your definition of indicatorFunction: if (abs(x[1]) > 0.25): return 1 #Phase1 else : - return 0 #Phase2 + return 2 #Phase2 + + +########### Options for material phases: ################################# +# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic" +######################################################################### +## Notation - Parameter input : +# 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 +###################################################################### +# --- Number of material phases +Phases=2 +#--- Define different material phases: +#- PHASE 1 +phase1_type="isotropic" +materialParameters_phase1 = [80, 80] + +#- PHASE 2 +phase2_type="isotropic" +materialParameters_phase2 = [80, 80] + + + +#--- define prestrain function for each phase +# (also works with non-constant values) +def prestrain_phase1(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def prestrain_phase2(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] diff --git a/materials/two_phase_material_3.py b/materials/two_phase_material_3.py index fea31cb1e031f7a2850503c19e0ff22a42d49b0a..f087f5d4c252b7bf5b75dcdbc1f0bbb91aef9e34 100644 --- a/materials/two_phase_material_3.py +++ b/materials/two_phase_material_3.py @@ -11,3 +11,34 @@ def f(x): return 1 #Phase1 else : return 0 #Phase2 + + +########### Options for material phases: ################################# +# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic" +######################################################################### +## Notation - Parameter input : +# 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 +###################################################################### +# --- Number of material phases +Phases=2 +#--- Define different material phases: +#- PHASE 1 +phase1_type="isotropic" +materialParameters_phase1 = [80, 80] + +#- PHASE 2 +phase2_type="isotropic" +materialParameters_phase2 = [80, 80] + + + +#--- define prestrain function for each phase +# (also works with non-constant values) +def prestrain_phase1(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] + +def prestrain_phase2(x): + return [[1, 0, 0], [0,1,0], [0,0,1]] diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc index af299a97ca38c14f7f4d8a35d41dde59d8c02057..9899928ebf8160562e8741fd0b7e37565869a10b 100644 --- a/src/Cell-Problem.cc +++ b/src/Cell-Problem.cc @@ -235,6 +235,7 @@ int main(int argc, char *argv[]) // --- get scale ratio double gamma = parameterSet.get<double>("gamma",1.0); + std::cout << "scale ratio (gamma) set to : " << gamma << std::endl; //------------------------------------------------------------------------------------------------ //--- Compute Correctors diff --git a/materials/material_backup.py b/src/deprecated_code/31-8-22_prestrainedMaterials/material_backup.py similarity index 100% rename from materials/material_backup.py rename to src/deprecated_code/31-8-22_prestrainedMaterials/material_backup.py