dune/microstructure/microstructure.hh

+#ifndef DUNE_MICROSTRUCTURE_HH
+#define DUNE_MICROSTRUCTURE_HH
+
+// add your classes here
+
+#endif // DUNE_MICROSTRUCTURE_HH

dune/microstructure/prestrainimp.hh

+#ifndef DUNE_MICROSTRUCTURE_PRESTRAINIMP_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAINIMP_HH
+
+// #include <dune/microstructure/microstructuredrodgrid.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+
+//template<class MicrostructuredRodGrid> reicht include
+
+using namespace Dune;
+
+
+class PrestrainImp {
+
+public:
+
+    static const int dim = 3;
+	static const int nCompo = 3;
+
+// 	using GridType_Ce = typename MicrostructuredRodGrid::CellGridType;
+    using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+    using GridView = CellGridType::LeafGridView;
+    using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+
+    using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+// 	using Domain = typename GridType_Ce::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+// 	using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+// 	using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+
+protected:
+
+	double p1, p2, theta;
+	double width; //cell geometry
+
+public:
+
+    PrestrainImp(double p1, double p2, double theta, double width) : p1(p1), p2(p2), theta(theta), width(width){}
+
+    Func2Tensor getPrestrain(unsigned int imp)
+    {
+    	using std::pow;
+	    using std::abs;
+	    using std::sqrt;
+	    using std::sin;
+	    using std::cos;
+
+    	if (imp==1)
+    	{
+
+            Func2Tensor B1_ = [this] (const Domain& x) {              //  ISOTROPIC PRESSURE
+                if (abs(x[0]) > (theta/2)  && x[2] > 0)
+                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                if (abs(x[0]) < (theta/2)  && x[2] < 0)
+                    return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+                else
+                    return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+
+            };
+            std::cout <<" Prestrain Type: 1 "<< std::endl;
+            return B1_;
+        }
+        else if (imp==2)
+        {
+            Func2Tensor B2_ = [this] (const Domain& x) {              // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE    
+
+                if (abs(x[0]) < (theta/2) && x[2] < 0 && x[2] > -(1.0/2.0) )
+                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                else
+                    return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+
+
+            };
+            std::cout <<" Prestrain Type: 2 "<< std::endl;
+            return B2_;
+
+
+        }
+
+
+
+
+            // TODO ANISOTROPIC PRESSURE
+
+//         else if (imp==2)
+//     	{
+//
+// 	    	Func2Tensor B2_ = [this] (const Domain& z)
+// 	    	{
+// 	    		auto pi = std::acos(-1.0);
+// 	      		double beta = pi/4.0 + pi/4.0*z[1]; //z[1]=x3
+// 	      		MatrixRT Id(0);
+// 	      		for (int i=0;i<nCompo;i++)
+// 	        		Id[i][i]=1.0;
+// 	      		MatrixRT pressure = Id;
+// 	      		pressure *= 1.0/6.0;
+// 	      		MatrixRT n_ot_n = {
+// 	      			{cos(beta)*cos(beta), 0.0, cos(beta)*sin(beta)},
+// 	      			{0.0,                           0.0, 0.0						  },
+// 	      			{cos(beta)*sin(beta), 0.0, sin(beta)*sin(beta)}
+// 	      		};
+// 	      		n_ot_n*=0.5;
+// 	      		pressure += n_ot_n;
+// 	      		pressure *= this->a;
+// 	      		return pressure;
+// 	    	};
+//     		return B2_;
+//     	}
+    }
+
+};
+
+
+#endif

inputs/cellsolver.parset

+# --- Parameter File as Input for 'Cell-Problem'
+#
+# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0
+# since otherwise these cant be read from other Files!
+# --------------------------------------------------------
+
+
+
+
+#path for logfile
+####outputPath = "../../outputs/output.txt"
+
+
+### Remove/Comment this when running via Python-Script:
+#outputPath = "../../outputs"
+
+
+
+#outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt"
+#outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs"
+
+
+#############################################
+#  Cell Domain
+#############################################
+# Domain 1: (-1/2, 1/2)^3  , Domain 2 : [0,1)^2 x (-1/2, 1/2)
+
+cellDomain = 1
+
+#############################################
+#  Grid parameters
+#############################################
+
+
+#######################################################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+########################################################################
+
+#numLevels =  1 3   # computes all levels from first to second entry
+#numLevels = 3 3     # computes all levels from first to second entry
+numLevels = 4 4    # computes all levels from first to second entry
+#numLevels = 1 6
+
+
+#Elements_Cell = 20 20 20	               # number elements in each direction (y1 y2 x3)
+#nElements_Cell = 30 30 30
+#nElements_Cell = 30 30 30
+#nElements_Cell = 50 50 50
+#nElements_Cell = 100 100 2
+#nElements_Cell = 100 100 100  // does not work
+#nElements_Cell = 10 10 10
+#nElements_Cell = 2 2 2
+#nElements_Cell = 4 4 4
+#nElements_Cell = 8 8 8
+#nElements_Cell = 16 16 16
+#nElements_Cell = 32 32 32
+#nElements_Cell = 64 64 64
+
+
+#gamma=50.0
+gamma=1.0
+#gamma=2.5
+
+#############################################
+#  Material parameters
+#############################################
+beta = 2.0    # ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case
+
+mu1=1.0
+lambda1=5.0
+
+
+
+rho1 = 1.0
+#alpha = 5.0    # ratio between prestrain parameters rho1 & rho2
+alpha = 2.0    # ratio between prestrain parameters rho1 & rho2
+
+
+theta = 0.125
+#theta = 0.25
+#theta = 0.3    # volume fraction   #default = 1.0/4.0
+#theta = 0.25   # volume fraction
+#theta = 0.75   # volume fraction
+
+
+
+
+
+
+material_prestrain_imp= "parametrized_Laminate"
+#material_prestrain_imp= "analytical_Example"
+#material_prestrain_imp ="isotropic_bilayer"
+#material_prestrain_imp= "circle_fiber"    #TEST
+
+
+# -- write Material and/or prestrain functions:
+
+#write_materialFunctions = true 
+#write_prestrainFunctions = true  # VTK norm of B ,
+
+#write_VTK = true
+
+
+
+# Prestrain Types:
+
+#1 Isotropic Pressure
+#    Func2Tensor B1_ = [this] (const Domain& x) {              // ISOTROPIC PRESSURE
+#                if (abs(x[0]) > (theta/2)  && x[2] > 0)
+#                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+#                if (abs(x[0]) < (theta/2)  && x[2] < 0)
+#                    return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+#                else
+#                    return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+#            };
+
+
+#############################################
+#  Assembly options
+#############################################
+
+set_IntegralZero = true
+#set_IntegralZero = false
+
+#arbitraryLocalIndex = 7
+#arbitraryElementNumber = 3
+
+#arbitraryLocalIndex = 0
+#arbitraryElementNumber = 0
+
+
+#############################################
+#  Solver Type
+
+Solvertype = 1
+
+Solver_verbosity = 0
+
+
+#write_corrector_phi1 = false
+#write_corrector_phi2 = false
+#write_corrector_phi3 = false
+#write_corrector_phi1 = true
+#write_corrector_phi2 = true
+#write_corrector_phi3 = true
+
+write_L2Error = true
+#write_IntegralMean = true
+
+
+#############################################
+#  Define Analytic Solutions

inputs/cellsolver_backup.parset

+# --- Parameter File as Input for 'Cell-Problem'
+#
+# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0
+# since otherwise these cant be read from other Files!
+# --------------------------------------------------------
+
+
+
+
+#path for logfile
+#outputPath = "../../outputs/output.txt"
+
+
+### Remove/Comment this when running via Python-Script:
+outputPath = "../../outputs"
+
+
+
+#outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt"
+#outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs"
+
+
+#############################################
+#  Cell Domain
+#############################################
+# Domain 1: (-1/2, 1/2)^3  , Domain 2 : [0,1)^2 x (-1/2, 1/2)
+
+cellDomain = 1
+
+#############################################
+#  Grid parameters
+#############################################
+
+
+#######################################################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+########################################################################
+
+#numLevels =  1 3   # computes all levels from first to second entry
+numLevels = 3 3     # computes all levels from first to second entry
+#numLevels = 1 6
+
+
+#Elements_Cell = 20 20 20	               # number elements in each direction (y1 y2 x3)
+#nElements_Cell = 30 30 30
+#nElements_Cell = 30 30 30
+#nElements_Cell = 50 50 50
+#nElements_Cell = 100 100 2
+#nElements_Cell = 100 100 100  // does not work
+#nElements_Cell = 10 10 10
+#nElements_Cell = 2 2 2
+#nElements_Cell = 4 4 4
+#nElements_Cell = 8 8 8
+#nElements_Cell = 16 16 16
+#nElements_Cell = 32 32 32
+#nElements_Cell = 64 64 64
+
+
+#gamma=50.0
+gamma=1.0
+#gamma=2.5
+
+#############################################
+#  Material parameters
+#############################################
+beta = 2.0    # ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case
+
+mu1=1.0
+lambda1=0.0
+#lambda1 = 20.0
+#lambda1 = 20.0
+#lambda1 = 5.0
+#mu1=1000.0
+
+
+rho1 = 1.0
+#alpha = 5.0    # ratio between prestrain parameters rho1 & rho2
+alpha = 2.0    # ratio between prestrain parameters rho1 & rho2
+
+
+theta = 0.125
+#theta = 0.25
+#theta = 0.3    # volume fraction   #default = 1.0/4.0
+#theta = 0.25   # volume fraction
+#theta = 0.75   # volume fraction
+
+####  material_implementation("analytical_Example") ?
+
+
+
+
+#material_prestrain_imp ="isotropic_bilayer"
+####  material_implementation("analytical_Example") ?
+
+material_prestrain_imp= "analytical_Example"
+#material_prestrain_imp ="isotropic_bilayer"
+
+
+#material_prestrain_imp= "circle_fiber"    #TEST
+write_prestrainFunctions = true  # VTK norm of B ,
+
+
+
+
+
+# Prestrain Types:
+
+#1 Isotropic Pressure
+#    Func2Tensor B1_ = [this] (const Domain& x) {              // ISOTROPIC PRESSURE
+#                if (abs(x[0]) > (theta/2)  && x[2] > 0)
+#                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+#                if (abs(x[0]) < (theta/2)  && x[2] < 0)
+#                    return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+#                else
+#                    return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+#            };
+
+
+#############################################
+#  Assembly options
+#############################################
+
+set_IntegralZero = true
+#set_IntegralZero = false
+
+#arbitraryLocalIndex = 7
+#arbitraryElementNumber = 3
+
+#arbitraryLocalIndex = 0
+#arbitraryElementNumber = 0
+
+
+#############################################
+#  Solver Type
+
+Solvertype = 1
+
+#write_corrector_phi1 = false
+#write_corrector_phi2 = false
+#write_corrector_phi3 = false
+#write_corrector_phi1 = true
+#write_corrector_phi2 = true
+#write_corrector_phi3 = true
+
+write_L2Error = true
+#write_IntegralMean = true
+
+
+#############################################
+#  Define Analytic Solutions

outputs/BMatrix.txt

+1 1 -0.324468879874192817
+1 2 -0.598155154475232376
+1 3 2.2946844494900733e-09

outputs/CellProblem-result-level1.vtu

+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="8" NumberOfPoints="27">
+      <PointData Vectors="Corrector phi_1 level1">
+        <DataArray type="Float32" Name="Corrector phi_1 level1" NumberOfComponents="3" format="ascii">
+          -4.92612e-17 5.12665e-16 -0.00520833 1.4585e-17 -2.06282e-16 -0.00520833 1.52502e-17 -9.02594e-17 -0.00520833 -1.6941e-18 -2.00537e-17 -0.00520833
+          -5.36618e-18 1.25837e-16 0.00520833 6.328e-18 -9.52493e-17 0.00520833 1.70486e-18 3.71619e-17 0.00520833 -1.5501e-18 -7.88106e-17 0.00520833
+          -4.92612e-17 5.12665e-16 -0.00520833 1.52502e-17 -9.02594e-17 -0.00520833 -5.36618e-18 1.25837e-16 0.00520833 1.70486e-18 3.71619e-17 0.00520833
+          -4.92612e-17 5.12665e-16 -0.00520833 1.4585e-17 -2.06282e-16 -0.00520833 -5.36618e-18 1.25837e-16 0.00520833 6.328e-18 -9.52493e-17 0.00520833
+          -4.92612e-17 5.12665e-16 -0.00520833 -5.36618e-18 1.25837e-16 0.00520833 8.50067e-18 3.84459e-17 -0.00520833 2.06099e-17 -1.12811e-16 -0.00520833
+          3.06229e-18 1.79722e-17 -0.00520833 -1.32859e-17 -1.17555e-16 -0.00520833 8.50067e-18 3.84459e-17 -0.00520833 3.06229e-18 1.79722e-17 -0.00520833
+          8.50067e-18 3.84459e-17 -0.00520833 2.06099e-17 -1.12811e-16 -0.00520833 8.50067e-18 3.84459e-17 -0.00520833
+        </DataArray>
+        <DataArray type="Float32" Name="Corrector phi_2 level1" NumberOfComponents="3" format="ascii">
+          -4.50589e-17 5.0976e-16 -0.00520833 1.13853e-17 -2.0196e-16 -0.00520833 1.94525e-17 -8.73542e-17 -0.00520833 -4.89382e-18 -2.43755e-17 -0.00520833
+          -5.01495e-18 1.25709e-16 0.00520833 5.86931e-18 -9.44243e-17 0.00520833 2.05608e-18 3.72899e-17 0.00520833 -2.00879e-18 -7.96356e-17 0.00520833
+          -4.50589e-17 5.0976e-16 -0.00520833 1.94525e-17 -8.73542e-17 -0.00520833 -5.01495e-18 1.25709e-16 0.00520833 2.05608e-18 3.72899e-17 0.00520833
+          -4.50589e-17 5.0976e-16 -0.00520833 1.13853e-17 -2.0196e-16 -0.00520833 -5.01495e-18 1.25709e-16 0.00520833 5.86931e-18 -9.44243e-17 0.00520833
+          -4.50589e-17 5.0976e-16 -0.00520833 -5.01495e-18 1.25709e-16 0.00520833 6.49215e-18 3.54358e-17 -0.00520833 2.18307e-17 -1.11602e-16 -0.00520833
+          1.05377e-18 2.09823e-17 -0.00520833 -1.2065e-17 -1.18764e-16 -0.00520833 6.49215e-18 3.54358e-17 -0.00520833 1.05377e-18 2.09823e-17 -0.00520833
+          6.49215e-18 3.54358e-17 -0.00520833 2.18307e-17 -1.11602e-16 -0.00520833 6.49215e-18 3.54358e-17 -0.00520833
+        </DataArray>
+        <DataArray type="Float32" Name="Corrector phi_3 level1" NumberOfComponents="3" format="ascii">
+          -2.67519e-19 4.0072e-18 -2.94019e-33 -3.45455e-19 -3.47827e-18 2.42632e-33 2.67519e-19 4.0072e-18 1.64678e-33 3.45455e-19 -3.47827e-18 -1.19935e-33
+          -2.80845e-19 1.75886e-18 -1.83096e-33 -5.65625e-19 -1.878e-18 1.60857e-33 2.80845e-19 1.75886e-18 7.48937e-34 5.65625e-19 -1.878e-18 -4.89971e-34
+          -2.67519e-19 4.0072e-18 -2.94019e-33 2.67519e-19 4.0072e-18 1.64678e-33 -2.80845e-19 1.75886e-18 -1.83096e-33 2.80845e-19 1.75886e-18 7.48937e-34
+          -2.67519e-19 4.0072e-18 -2.94019e-33 -3.45455e-19 -3.47827e-18 2.42632e-33 -2.80845e-19 1.75886e-18 -1.83096e-33 -5.65625e-19 -1.878e-18 1.60857e-33
+          -2.67519e-19 4.0072e-18 -2.94019e-33 -2.80845e-19 1.75886e-18 -1.83096e-33 2.05575e-19 1.16377e-18 -1.54544e-33 -1.79765e-18 -1.45444e-18 1.35312e-33
+          -2.05575e-19 1.16377e-18 4.59892e-34 1.79765e-18 -1.45444e-18 -2.74301e-34 2.05575e-19 1.16377e-18 -1.54544e-33 -2.05575e-19 1.16377e-18 4.59892e-34
+          2.05575e-19 1.16377e-18 -1.54544e-33 -1.79765e-18 -1.45444e-18 1.35312e-33 2.05575e-19 1.16377e-18 -1.54544e-33
+        </DataArray>
+      </PointData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          -0.5 -0.5 -0.5 0 -0.5 -0.5 -0.5 0 -0.5 0 0 -0.5
+          -0.5 -0.5 0 0 -0.5 0 -0.5 0 0 0 0 0
+          0.5 -0.5 -0.5 0.5 0 -0.5 0.5 -0.5 0 0.5 0 0
+          -0.5 0.5 -0.5 0 0.5 -0.5 -0.5 0.5 0 0 0.5 0
+          0.5 0.5 -0.5 0.5 0.5 0 -0.5 -0.5 0.5 0 -0.5 0.5
+          -0.5 0 0.5 0 0 0.5 0.5 -0.5 0.5 0.5 0 0.5
+          -0.5 0.5 0.5 0 0.5 0.5 0.5 0.5 0.5
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 4 5 7 6 1 8 9 3
+          5 10 11 7 2 3 13 12 6 7 15 14
+          3 9 16 13 7 11 17 15 4 5 7 6
+          18 19 21 20 5 10 11 7 19 22 23 21
+          6 7 15 14 20 21 25 24 7 11 17 15
+          21 23 26 25
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          8 16 24 32 40 48 56 64
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          12 12 12 12 12 12 12 12
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>