Matlab-Programs/symMinimization.m

+
+function [G, angle, Type, kappa] = symMinimization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath) %(Q_hom,B_eff)
+
+
+syms v1 v2 q1 q2 q3 q12 b1 b2 b3
+
+
+% -------- Options ----------
+% % print_Input = false;  %effective quantities
+% print_Input = true;
+% % print_statPoint = false;
+% print_statPoint = true;
+% print_Minimizer = true;
+% % make_FunctionPlot = false;
+% make_FunctionPlot = true;
+
+print_Uniqueness = false;         
+check_StationaryPoints = false;   % could be added to Input-Parameters..
+compare_with_Classification = false;  %maybe added later --- 
+
+
+%fprintf('Functions to be minimized:')
+f_plus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2-2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2)...
+                                        + q12*(v1^2*v2^2-b2*v1^2-b1*v2^2+b1*b2);
+f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2)...
+                                        + q12*(v1^2*v2^2+b2*v1^2+b1*v2^2+b1*b2);
+
+
+% ---- Fix parameters
+
+
+% Epsilon used:
+epsilon = 1e-8;
+
+
+if ~exist('InputPath','var')
+     % third parameter does not exist, so default it to something
+      absPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs";
+ end
+
+
+% 1. Import effective quantities from CellSolver-Code:
+%read as sparse Matrix...
+try %absolutePath
+    Qmat = spconvert(load(absPath + '' + "/QMatrix.txt"));
+    Bmat = spconvert(load(absPath + '' + "/BMatrix.txt"));
+%     fprintf('Use absolute Path')
+catch ME  % use relativePath
+    Qmat = spconvert(load('../outputs/QMatrix.txt'));
+    Bmat = spconvert(load('../outputs/BMatrix.txt'));
+%     fprintf('Use relative Path')
+end
+%convert to full matrix... 
+Qmat = full(Qmat);
+Bmat = full(Bmat);
+
+
+
+
+% --- TODO CHECK: assert if Q is orthotropic ???  check ifq13=q31=q23=q32= 0 ?
+
+
+if print_Input
+    fprintf('effective quadratic form:')
+    Qmat
+    fprintf('effective prestrain')
+    Bmat
+    
+    % check if Q is (close to..) symmetric
+    % könnte Anti-symmetric part berechnen und schauen dass dieser klein?
+    % Test: issymmetric(Qmat) does not work for float matrices?
+    % symmetric part 0.5*(Qmat+Qmat')
+    % anti-symmetric part 0.5*(Qmat-Qmat')
+    if norm(0.5*(Qmat-Qmat'),'fro') < 1e-8
+        fprintf('Qmat (close to) symmetric \n')
+        norm(0.5*(Qmat-Qmat'),'fro') % TEST
+    else
+        fprintf('Qmat not symmetric \n')
+    end
+    
+    
+    % Check if B_eff is diagonal this is equivalent to b3 == 0
+    if abs(Bmat(3)) < 1e-8
+        fprintf('B_eff is diagonal (b3 == 0) \n')
+    else
+        fprintf('B_eff is  NOT diagonal (b3 != 0)  \n')
+    end
+end
+
+
+
+
+% CAST VALUES TO SYM FIRST? This is done anyway..
+% % Substitute effective quantitites      
+f_plus = subs(f_plus,{q1, q2, q3, q12, b1, b2, b3}, {Qmat(1,1), Qmat(2,2), Qmat(3,3), Qmat(1,2), ...
+                                                     Bmat(1), Bmat(2), Bmat(3)});
+
+f_minus = subs(f_minus,{q1, q2, q3, q12, b1, b2, b3}, {Qmat(1,1), Qmat(2,2), Qmat(3,3), Qmat(1,2), ...
+                                                     Bmat(1), Bmat(2), Bmat(3)});
+
+
+% Compute the Gradients 
+df_plusx = diff(f_plus,v1);
+df_plusy = diff(f_plus,v2);
+df_minusx = diff(f_minus,v1);
+df_minusy = diff(f_minus,v2);
+
+% Setup Equations Grad(f) = 0 
+eq1 = df_plusx == 0;
+eq2 = df_plusy == 0;
+eqns_plus = [eq1, eq2];
+eq3 = df_minusx == 0;
+eq4 = df_minusy == 0;
+eqns_minus = [eq3, eq4];
+
+% -------  Symbolically Solve Equations: 
+% More robust (works even for values b_3 ~ 1e-08 ): 
+S_plus = solve(eqns_plus,v1,v2,'MaxDegree' , 5);  
+S_minus = solve(eqns_minus,v1,v2,'MaxDegree' , 5);
+
+A_plus = S_plus.v1;
+B_plus = S_plus.v2;
+A_minus = S_minus.v1;
+B_minus = S_minus.v2;
+
+if check_StationaryPoints
+    %---------- TEST if Grad(f) = 0 ---------------------
+    fprintf('Testing equation grad(f) = 0  with stationary points')
+    for i = 1:size(A_plus,1)
+        fprintf('Testing %d.point (f_plus): ',i )
+        [ double(subs(subs(df_plusx,v1,A_plus(i)),v2,B_plus(i))), double(subs(subs(df_plusy,v1,A_plus(i)),v2,B_plus(i))) ]
+    end
+    for i = 1:size(A_minus,1)
+        fprintf('Testing %d.point (f_minus): ',i )
+        [double(subs(subs(df_minusx,v1,A_minus(i)),v2,B_minus(i))), double(subs(subs(df_minusy,v1,A_minus(i)),v2,B_minus(i)))]
+    end
+    % ------------------------------------
+end
+
+% --- Extract only Real-Solutions
+% fprintf('real stationary points of f_plus:')
+tmp1 = A_plus(imag(double(A_plus))==0 & imag(double(B_plus)) == 0);
+tmp2 = B_plus(imag(double(A_plus))==0 & imag(double(B_plus)) == 0);
+A_plus = tmp1;
+B_plus = tmp2;
+SP_plus = [A_plus,B_plus];
+
+% fprintf('real stationary points of f_minus:')
+tmp1 = A_minus(imag(double(A_minus))==0 & imag(double(B_minus)) == 0);
+tmp2 = B_minus(imag(double(A_minus))==0 & imag(double(B_minus)) == 0);
+A_minus = tmp1;
+B_minus = tmp2;
+SP_minus = [A_minus,B_minus];
+
+% TODO one should use f_plus.subs(A_plus..) to compute function value symbolically?
+% in the end only the stationaryPoints are used.. Ok to compare function values numerically
+
+% Determine global Minimizer from stationary points:
+% fprintf('function values at stationary points (f_plus):')
+T_plus = arrayfun(@(v1,v2) double(f_plus(v1,v2,q1,q2,q3,q12,b1,b2,b3)),A_plus,B_plus,'UniformOutput', false);
+T_plus = cell2mat(T_plus);
+
+%Test: use Substitution
+% subs(f_plus,{v1, v2}, {A_plus,B_plus})
+
+% fprintf('function values at stationary points (f_minus):')
+T_minus = arrayfun(@(v1,v2) double(f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3)),A_minus,B_minus,'UniformOutput', false);
+T_minus = cell2mat(T_minus);
+
+%Test: use Substitution
+% T_minus = subs(f_minus,{v1, v2}, {A_minus,B_minus})
+% double(T_minus)
+
+if print_statPoint
+    fprintf('real stationary points of f_plus: ')
+%     SP_Plus     %alternative: output as symbolic (can be unwieldy)
+    double(SP_plus)
+    fprintf('real stationary points of f_minus:')
+%     SP_Minus     %alternative: output as symbolic (can be unwieldy)
+    double(SP_minus)
+    fprintf('function values at stationary points (f_plus):')
+    T_plus
+    fprintf('function values at stationary points (f_minus):')
+    T_minus
+end
+
+% --- Find Stationary Point(s) with minimum Value 
+[Min_plus,MinIdx_plus] = min(T_plus, [], 'all', 'linear');    %find one min...
+[Min_minus,MinIdx_minus] = min(T_minus, [], 'all', 'linear');
+% [Min_minus,MinIdx_minus] = min(T_minus)                     % works with symbolic too
+
+% Compare Minimizers of f_plus & f_minus
+[globalMinimizerValue,GlobalIdx] = min([Min_plus,Min_minus]); 
+
+if GlobalIdx == 1     %Min_plus   % i.e. Global Minimizer is given by f_plus
+    GlobalMinimizer = SP_plus(MinIdx_plus,:);
+    sign = 1.0;
+elseif GlobalIdx == 2  %Min_minus % i.e. Global Minimizer is given by f_minus
+    GlobalMinimizer = SP_minus(MinIdx_minus,:);
+    sign = -1.0; 
+end
+
+% ------ Check if there are more SP with the same value...
+ MinIndices_minus = find(T_minus(:) == globalMinimizerValue);    % Find indices of All Minima
+ MinIndices_plus  = find(T_plus(:) == globalMinimizerValue);    % Find indices of All Minima                   
+ numMinSP_minus = size(MinIndices_minus,1);   % One of these is always >= 2 due to the structure of the roots..
+ numMinSP_plus  = size(MinIndices_plus,1);
+    
+%     AllMinSP_minus = SP_minus(MinIndices_minus,:)
+%     AllMinSP_minus = double(SP_minus(MinIndices_minus,:))
+%     AllMin = T_minus(MinIndices)  %bereits klar dass diese selben funktionswert haben..
+      
+    Minimizer = sign*(GlobalMinimizer'*GlobalMinimizer);  % global minimizing Matrix G*
+    MinimizerCount = 1;
+    % different Stationary Points might correspond to the same minimizing
+    % Matrix G*... check this: 
+    
+    % Compare only with other StationaryPoints/Minimizers
+    % remove Index of Minimizer
+    if GlobalIdx == 1     
+        MinIndices_plus = MinIndices_plus(MinIndices_plus~=MinIdx_plus);
+    elseif GlobalIdx == 2  
+        MinIndices_minus = MinIndices_minus(MinIndices_minus~=MinIdx_minus);
+    end
+    
+    MinIndices = cat(1,MinIndices_plus,MinIndices_minus);  %[Minimizers-Indices f_plus, Minimizer-Indices f_minus]
+    
+    for i = 1:(numMinSP_minus+numMinSP_plus-1)  % -1: dont count Minimizer itself.. 
+        idx = MinIndices(i);   
+        
+        if i > numMinSP_plus
+            SP = SP_minus(idx,:);
+        else
+            SP = SP_plus(idx,:); 
+        end
+        
+%         SP_value = T_minus(idx) % not needed? 
+        Matrix = sign*(SP'*SP);
+
+        if norm(double(Matrix-Minimizer),'fro') < 1e-8   %check is this sufficient here?
+%             fprintf('both StationaryPoints correspond to the same(Matrix-)Minimizer')
+        else
+%             fprintf('StationaryPoint corresponds to a different (Matrix-)Minimizer')
+            MinimizerCount = MinimizerCount + 1;
+        end
+    end
+% ----------------------------------------------------------------------------------------------------------------
+
+
+
+% Output Uniqueness of Minimizers:
+if print_Uniqueness
+    if MinimizerCount == 1
+        fprintf('Unique Minimzier')
+    elseif MinimizerCount == 2
+        fprintf('Two Minimziers')
+    else
+        fprintf('1-Parameter family of Minimziers')
+    end
+end
+
+
+
+
+
+% --- determine the angle of the Minimizer
+% a1 = Minimizer(1,1)
+% a2 = Minimizer(2,2)
+a1 = double(Minimizer(1,1));
+a2 = double(Minimizer(2,2));
+
+% compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
+e = [sqrt((a1/(a1+a2))), sqrt((a2/(a1+a2)))];    % always positive under sqrt here .. basically takes absolute value here
+angle = atan2(e(2), e(1)); 
+
+% compute curvature kappa
+kappa = (a1 + a2);
+
+% % CHeck off diagonal entries:
+% sqrt(a1*a2);
+% double(Minimizer);
+
+G = double(Minimizer);
+
+
+
+% --- "Classification" / Determine the TYPE of Minimizer by using 
+% the number of solutions (Uniqueness?)
+% the angle (axial- or non-axial Minimizer)
+
+% (Alternative compute det[GlobalMinimizer' e1']  where e1 = [1 0] ?)
+
+% Check Uniqueness -- Options: unique/twoMinimizers/1-ParameterFamily
+if MinimizerCount == 1
+%     fprintf('Unique Minimzier')
+    % Check if Minimizer is axial or non-axial:
+    if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
+     Type = 3;
+    else % non-axial Minimizer   
+     Type = 1;
+    end
+elseif MinimizerCount == 2
+%     fprintf('Two Minimziers')
+    % Check if Minimizer is axial or non-axial:
+    if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
+        Type = 3;
+    else % non-axial Minimizer   
+        fprintf('ERROR: Two non-axial Minimizers cannot happen!')
+    end
+else
+%     fprintf('1-Parameter family of Minimziers')
+    % Check if Minimizer is axial or non-axial:
+    if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
+%         fprintf('ERROR: axial Minimizers cannot happen for 1-Parameter Family!')
+    else % non-axial Minimizer   
+        Type = 2;
+    end
+end
+% ------------------------------------------------------------------------------------------------------
+
+
+
+
+if print_Output
+    fprintf(' --------- Output symMinimization --------')
+    fprintf('Global Minimizer v: (%d,%d) \n', GlobalMinimizer(1),GlobalMinimizer(2) )
+    fprintf('Global Minimizer Value f(v): %d \n', sym(globalMinimizerValue) ) %cast to symbolic
+    %     fprintf('Global Minimizer Value : %d', globalMinimizerValue )
+    
+    fprintf('Global Minimizer G: \n' )
+    G 
+    fprintf("Angle = %d \n", angle)
+    fprintf("Curvature = %d \n", kappa)
+    fprintf("Type = %i \n", Type)
+    fprintf(' --------- -------------------- --------')
+end
+
+
+
+if make_FunctionPlot  
+    fsurf(@(x,y) f_plus(x,y,q1,q2,q3,q12,b1,b2,b3)) % Plot functions
+    hold on
+    plot3(double(A_plus),double(B_plus),T_plus,'g*')
+    %Plot GlobalMinimizer:
+    hold on
+    plot3(double(GlobalMinimizer(1)),double(GlobalMinimizer(2)),globalMinimizerValue, 'o', 'Color','c')
+    % view(90,0)
+    % view(2)
+
+    figure
+    fsurf(@(x,y) f_minus(x,y,q1,q2,q3,q12,b1,b2,b3))
+    hold on
+    plot3(double(A_minus),double(B_minus),T_minus,'g*')
+    hold on
+    plot3(double(GlobalMinimizer(1)), double(GlobalMinimizer(2)),globalMinimizerValue, 'o', 'Color','c')
+end
+
+
+return 
+
+
+% Write symbolic solution to txt-File in Latex format
+% fileID = fopen('txt.txt','w');
+% fprintf(fileID,'%s' , latex(S_plus.v1));
+% fclose(fileID);
+
+
+
+
+
+
+

README

+Preparing the Sources
+=========================
+
+Additional to the software mentioned in README you'll need the
+following programs installed on your system:
+
+  cmake >= 3.1
+
+Getting started
+---------------
+
+If these preliminaries are met, you should run
+
+  dunecontrol all
+
+which will find all installed dune modules as well as all dune modules
+(not installed) which sources reside in a subdirectory of the current
+directory. Note that if dune is not installed properly you will either
+have to add the directory where the dunecontrol script resides (probably
+./dune-common/bin) to your path or specify the relative path of the script.
+
+Most probably you'll have to provide additional information to dunecontrol
+(e. g. compilers, configure options) and/or make options.
+
+The most convenient way is to use options files in this case. The files
+define four variables:
+
+CMAKE_FLAGS      flags passed to cmake (during configure)
+
+An example options file might look like this:
+
+#use this options to configure and make if no other options are given
+CMAKE_FLAGS=" \
+-DCMAKE_CXX_COMPILER=g++-5 \
+-DCMAKE_CXX_FLAGS='-Wall -pedantic' \
+-DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-5 and set compiler flags
+
+If you save this information into example.opts you can pass the opts file to
+dunecontrol via the --opts option, e. g.
+
+  dunecontrol --opts=example.opts all
+
+More info
+---------
+
+See
+
+     dunecontrol --help
+
+for further options.
+
+
+The full build system is described in the dune-common/doc/buildsystem (Git version) or under share/doc/dune-common/buildsystem if you installed DUNE!

config.h.cmake

+/* begin dune-microstructure
+   put the definitions for config.h specific to
+   your project here. Everything above will be
+   overwritten
+*/
+
+/* begin private */
+/* Name of package */
+#define PACKAGE "@DUNE_MOD_NAME@"
+
+/* Define to the address where bug reports for this package should be sent. */
+#define PACKAGE_BUGREPORT "@DUNE_MAINTAINER@"
+
+/* Define to the full name of this package. */
+#define PACKAGE_NAME "@DUNE_MOD_NAME@"
+
+/* Define to the full name and version of this package. */
+#define PACKAGE_STRING "@DUNE_MOD_NAME@ @DUNE_MOD_VERSION@"
+
+/* Define to the one symbol short name of this package. */
+#define PACKAGE_TARNAME "@DUNE_MOD_NAME@"
+
+/* Define to the home page for this package. */
+#define PACKAGE_URL "@DUNE_MOD_URL@"
+
+/* Define to the version of this package. */
+#define PACKAGE_VERSION "@DUNE_MOD_VERSION@"
+
+/* end private */
+
+/* Define to the version of dune-microstructure */
+#define DUNE_MICROSTRUCTURE_VERSION "@DUNE_MICROSTRUCTURE_VERSION@"
+
+/* Define to the major version of dune-microstructure */
+#define DUNE_MICROSTRUCTURE_VERSION_MAJOR @DUNE_MICROSTRUCTURE_VERSION_MAJOR@
+
+/* Define to the minor version of dune-microstructure */
+#define DUNE_MICROSTRUCTURE_VERSION_MINOR @DUNE_MICROSTRUCTURE_VERSION_MINOR@
+
+/* Define to the revision of dune-microstructure */
+#define DUNE_MICROSTRUCTURE_VERSION_REVISION @DUNE_MICROSTRUCTURE_VERSION_REVISION@
+
+/* end dune-microstructure
+   Everything below here will be overwritten
+*/

doc/CMakeLists.txt

+add_subdirectory("doxygen")

doc/doxygen/CMakeLists.txt

+# shortcut for creating the Doxyfile.in and Doxyfile
+add_doxygen_target()

doc/doxygen/Doxylocal

+# This file contains local changes to the doxygen configuration
+# please use '+=' to add files/directories to the lists
+
+# The INPUT tag can be used to specify the files and/or directories that contain
+# documented source files. You may enter file names like "myfile.cpp" or
+# directories like "/usr/src/myproject". Separate the files or directories
+# with spaces.
+
+INPUT                 += @top_srcdir@/dune/
+# see e.g. dune-grid for the examples of mainpage and modules
+# INPUT                 += @srcdir@/mainpage \
+#                          @srcdir@/modules
+
+# The EXCLUDE tag can be used to specify files and/or directories that should
+# be excluded from the INPUT source files. This way you can easily exclude a
+# subdirectory from a directory tree whose root is specified with the INPUT tag.
+
+# EXCLUDE               += @top_srcdir@/dune/microstructure/test
+
+# The EXAMPLE_PATH tag can be used to specify one or more files or
+# directories that contain example code fragments that are included (see
+# the \include command).
+
+# EXAMPLE_PATH          += @top_srcdir@/src
+
+# The IMAGE_PATH tag can be used to specify one or more files or
+# directories that contain image that are included in the documentation (see
+# the \image command).
+
+# IMAGE_PATH            += @top_srcdir@/dune/microstructure/pics

dune-microstructure.pc.in

+prefix=@prefix@
+exec_prefix=@exec_prefix@
+libdir=@libdir@
+includedir=@includedir@
+CXX=@CXX@
+CC=@CC@
+DEPENDENCIES=@REQUIRES@
+
+Name: @PACKAGE_NAME@
+Version: @VERSION@
+Description: dune-microstructure module
+URL: http://dune-project.org/
+Requires: dune-common dune-istl dune-geometry dune-localfunctions dune-uggrid dune-grid dune-typetree dune-functions dune-solvers
+Libs: -L${libdir}
+Cflags: -I${includedir}

dune.module

+################################
+# Dune module information file #
+################################
+
+# Name of the module
+Module: dune-microstructure
+Version: 1.0
+Maintainer: klaus.boehnlein@tu-dresden.de
+# Required build dependencies
+Depends: dune-common dune-istl dune-geometry dune-localfunctions dune-uggrid dune-grid dune-typetree dune-functions 
+# Optional build dependencies
+#Suggests:

dune/CMakeLists.txt

+add_subdirectory(microstructure)

dune/microstructure/CMakeLists.txt

+#install headers
+install(FILES microstructure.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/microstructure)

dune/microstructure/matrix_operations.hh

+#ifndef DUNE_MICROSTRUCTURE_MATRIXOPERATIONS_HH
+#define DUNE_MICROSTRUCTURE_MATRIXOPERATIONS_HH
+
+namespace MatrixOperations {
+
+	using namespace Dune;
+
+	using MatrixRT = FieldMatrix< double, 3, 3>;
+	using VectorRT = FieldVector< double, 3>;	
+
+	using std::sin;
+	using std::cos;
+
+	static MatrixRT Id (){ 
+		MatrixRT Id(0);
+	    for (int i=0;i<3;i++)
+	        Id[i][i]=1.0;
+	    return Id;
+	}
+	
+	
+    static double norm(VectorRT v){
+		return sqrt(pow(v[0],2) + pow(v[1],2) + pow(v[2],2));
+	}
+
+	static double norm(MatrixRT M){
+		return sqrt(
+			  pow(M[0][0],2) + pow(M[0][1],2) + pow(M[0][2],2)
+			+ pow(M[1][0],2) + pow(M[1][1],2) + pow(M[1][2],2)
+			+ pow(M[2][0],2) + pow(M[2][1],2) + pow(M[2][2],2));
+	}
+	
+
+	static MatrixRT sym (MatrixRT M) { // 1/2 (M^T + M)
+    	MatrixRT ret(0);
+    	for (int i = 0; i< 3; i++)
+    		for (int j = 0; j< 3; j++)
+    			ret[i][j] = 0.5 * (M[i][j] + M[j][i]);
+    	return ret;
+	}
+
+	static VectorRT crossProduct (VectorRT v, VectorRT w) { // v otimes w
+    	return {v[1]*w[2] - v[2]*w[1], -1*(v[0]*w[2] - v[2]*w[0]), v[0]*w[1] - v[1]*w[0]};
+	}
+
+	static MatrixRT rankoneTensorproduct (VectorRT v, VectorRT w) { // v otimes w
+    	return 
+    	{{v[0]*w[0], v[1]*w[0], v[2]*w[0]}, 
+	     {v[0]*w[1], v[1]*w[1], v[2]*w[1]}, 
+	     {v[0]*w[2], v[1]*w[2], v[2]*w[2]}};
+	}
+
+	static MatrixRT nematicLiquidCrystal (double p, VectorRT n){ //B = 1/6*p*Id + 1/2*p*(n otimes n)
+		MatrixRT B(0);
+	    for (int i=0;i<3;i++)
+	        B[i][i]=p/6.0;
+	    MatrixRT n_ot_n = rankoneTensorproduct(n,n);
+	    n_ot_n*=p/2.0;
+	    B += n_ot_n;
+	    return B;
+	}
+
+	static MatrixRT biotStrainApprox (VectorRT U, VectorRT k, VectorRT e_cs){ //E_h = (U +  k x e_cs, 0, 0)
+	    VectorRT k_x_ecs = crossProduct(k, e_cs);
+	    VectorRT U_plus_k_x_ecs = U + k_x_ecs;
+	    VectorRT e_1 = {1, 0, 0};
+	    return rankoneTensorproduct(U_plus_k_x_ecs, e_1);
+	}
+
+	static MatrixRT crossSectionDirectionScaling(double w, MatrixRT M){
+		return {M[0], w*M[1], w*M[2]};
+	}
+
+	static double trace (MatrixRT M){ 
+	    return M[0][0]+ M[1][1] + M[2][2];
+	}
+
+	static double scalarProduct (MatrixRT M1, MatrixRT M2){ 
+		double sum = 0;
+		for (int i=0; i<3; i++)
+		    for (int j=0; j<3; j++)
+		    	sum += M1[i][j] * M2[i][j];
+	    return sum;
+	}
+
+// 	static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){    // ?? Check with Robert
+// 		E1= sym(E1);
+// 		E2 = sym(E2);
+// 		double t1 = scalarProduct(E1,E2);
+// 		t1 *= 2* mu;
+// 		double t2 = trace(E1)*trace(E2); 
+// 		t2 *= lambda;
+// 		return t1 + t2;
+// 	}
+	
+    static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2)  // CHANGED
+    {  
+        auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+        return scalarProduct(t1,E2);
+
+	}
+
+    // --- Generalization: Define Quadratic QuadraticForm
+    
+    static double QuadraticForm(const double mu, const double lambda, const MatrixRT M){
+        
+        auto tmp1 = sym(M);
+        double tmp2 = norm(tmp1);
+        return lambda*std::pow(trace(M),2) + 2*mu*pow( tmp2 ,2);
+//         return lambda*std::pow(trace(M),2) + 2*mu*pow( norm( sym(M) ),2);
+    }
+
+    
+	static double generalizedDensity(const double mu, const double lambda, MatrixRT F, MatrixRT G){
+     /// Write this whole File as a Class that uses lambda,mu as members ? 
+        
+     // Define L via Polarization-Identity from QuadratifForm
+     // <LF,G> := (1/2)*(Q(F+G) - Q(F) - Q(G) ) 
+        return (1.0/2.0)*(QuadraticForm(mu,lambda,F+G) - QuadraticForm(mu,lambda,F) - QuadraticForm(mu,lambda,G) );
+    }
+
+	static MatrixRT matrixSqrt(MatrixRT M){
+		std::cout << "matrixSqrt not implemented!!!" << std::endl;//implement this
+		return M;
+	}
+
+	static double lameMu(double E, double nu){
+		return 0.5 * 1.0/(1.0 + nu) * E;
+	}
+
+	static double lameLambda(double E, double nu){
+		return nu/(1.0-2.0*nu) * 1.0/(1.0+nu) * E;
+	}
+
+	static bool isInRotatedPlane(double phi, double x1, double x2){
+		return cos(phi)*x1 + sin(phi)*x2 > 0;
+	}
+
+
+
+	/*
+	template<double phi>
+	static bool isInRotatedPlane(double x1, double x2){
+		return cos(phi)*x1 + sin(phi)*x2 > 0;
+	}*/
+
+}
+
+
+
+// OPTIONAL H1-Seminorm ... 
+
+/*
+template<class VectorCT>
+double semi_H1_vectorScalarProduct(VectorCT& coeffVector1, VectorCT& coeffVector2)
+{
+
+double ret(0);
+auto localView = basis_.localView();
+
+for (const auto& e : elements(basis_.gridView()))
+  	{
+	    localView.bind(e);
+
+	    auto geometry = e.geometry();
+
+	    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+	    const auto nSf = localFiniteElement.localBasis().size();
+
+	    int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+	    const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
+
+	   	for (const auto& quadPoint : quad)
+    	{
+    		double elementScp(0);
+
+    		auto pos = quadPoint.position();
+    		const auto jacobian = geometry.jacobianInverseTransposed(pos);
+    		const double integrationElement = geometry.integrationElement(pos);
+
+    		std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+    		localFiniteElement.localBasis().evaluateJacobian(pos, referenceGradients);
+
+    		std::vector< VectorRT > gradients(referenceGradients.size());
+    		for (size_t i=0; i< gradients.size(); i++)
+      		jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+    		for (size_t i=0; i < nSf; i++)
+    			for (size_t j=0; j < nSf; j++)
+    			for (size_t k=0; k < dimworld; k++)
+    				for (size_t l=0; l < dimworld; l++)
+          		{
+          			MatrixRT defGradient1(0);
+          			defGradient1[k] = gradients[i];
+
+          			MatrixRT defGradient2(0);
+          			defGradient2[l] = gradients[j];
+
+		          	size_t localIdx1 = localView.tree().child(k).localIndex(i);
+		          	size_t globalIdx1 = localView.index(localIdx1);
+
+		          	size_t localIdx2 = localView.tree().child(l).localIndex(j);
+		          	size_t globalIdx2 = localView.index(localIdx2);
+
+        			double scp(0);
+        			for (int ii=0; ii<dimworld; ii++)
+            		for (int jj=0; jj<dimworld; jj++)
+        						scp += coeffVector1[globalIdx1] * defGradient1[ii][jj] * coeffVector2[globalIdx2] * defGradient2[ii][jj];
+
+        				elementScp += scp;
+        			}
+
+            ret += elementScp * quadPoint.weight() * integrationElement;
+      } //qp
+  	} //e
+	  return ret;
+}*/
+
+
+
+
+
+
+/*
+	class BiotStrainApprox
+	{
+
+		// E_h = h^{-1} (R^T F_h - Id)
+		//
+		//E_h = (U +  K e_cs, 0, 0)
+		//
+
+		Func2Tensor E_a = [] (const Domain& z) //extra Klasse
+	  	{
+			return biotStrainApprox({1, 0, 0}, {0, 0, 0}, {0, z[1], z[2]}); //MatrixRT{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; 
+	  	};
+
+	  	// twist in x2-x3 plane E_K1 = sym (e1 x (0,x2,x3)^T ot e1)
+	    Func2Tensor E_K1 = [] (const Domain& z) 
+	    {
+	    	return biotStrainApprox({0, 0, 0}, {1, 0, 0}, {0, z[1], z[2]}); //MatrixRT{{0.0,-0.5*z[2], 0.5*z[1]}, {-0.5*z[2], 0.0 , 0.0 }, {0.5*z[1],0.0,0.0}}; 
+	  	};
+
+	    //  bend strain in x1-x2 plane	E_K2 = sym (e2 x (0,x2,x3)^T ot e1)
+	  	Func2Tensor E_K2 = [] (const Domain& z) 
+	  	{
+	  		return biotStrainApprox({0, 0, 0}, {0, 1, 0}, {0, z[1], z[2]}); //MatrixRT{{z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0}}; 
+	  	};
+
+	  	// bend strain in x1-x3 plane	E_K3 = sym (e3 x (0,x2,x3)^T ot e1)
+	  	Func2Tensor E_K3 = [] (const Domain& z) 
+	  	{
+	  		return biotStrainApprox({0, 0, 0}, {0, 0, 1}, {0, z[1], z[2]}); //MatrixRT{{-z[1], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; 
+	  	};
+	};
+*/
+
+#endif

dune/microstructure/microstructure.hh

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

dune/microstructure/prestrain_material_geometry.hh

+#ifndef DUNE_MICROSTRUCTURE_PRESTRAIN_MATERIAL_GEOMETRY_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAIN_MATERIAL_GEOMETRY_HH
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/microstructure/matrix_operations.hh>
+
+using namespace Dune;
+using namespace MatrixOperations;
+using std::pow;
+using std::abs;
+using std::sqrt;
+using std::sin;
+using std::cos;
+
+
+
+
+
+template <int dim>
+class IsotropicMaterialImp 
+{
+
+public:
+	static const int dimWorld = dim;
+    using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+    using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+    using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+	using FuncScalar = std::function< double(const Domain&) >;
+
+
+    FuncScalar getMu(ParameterTree parameters){
+
+    	std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
+    	//std::array<std::string,5> imps({"homogen_poisson" "bilayer_poisson" "chess_poisson" "3Dchess_poisson" "vertical_bilayer_poisson"});
+        
+    	double phi = M_PI*parameters.get<double>("material_angle", 0.0)/180; //TODO 
+    	double theta = parameters.get<double>("theta",1.0/4.0);                        //TODO alle parameter hier laden?
+        
+        double width = parameters.get<double>("width", 1.0);
+        double height = parameters.get<double>("height", 1.0);
+
+
+    	if (imp == "homogeneous"){    
+		    double mu1     = parameters.get<double>("mu1", 10);
+
+		    auto muTerm = [mu1] (const Domain& z) {return mu1;};
+		    
+			return muTerm;
+		}
+		else if (imp == "isotropic_bilayer")
+    	{
+            double mu1     = parameters.get<double>("mu1",10.0);
+            double beta = parameters.get<double>("beta",2.0); 
+		    double mu2 = beta*mu1;
+    
+ 		   auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+                    if (z[2]>0)                                                   
+                        return mu1;
+                    else
+                        return mu2;
+                    };
+		    
+           return muTerm;
+		}
+        else if (imp == "analytical_Example"){	
+		    double mu1     = parameters.get<double>("mu1",10.0);
+            double beta = parameters.get<double>("beta",2.0); 
+		    double mu2 = beta*mu1;
+		    
+		   auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+//                     std::cout << "Analytical-MU is used" << std::endl;
+                    if (abs(z[0]) >= (theta/2.0))                                                   
+                        return mu1;
+                    else
+                        return mu2;
+                    };
+		    
+           return muTerm;
+		}
+		else if (imp == "circle_fiber"){
+		    double mu1     = parameters.get<double>("mu1",10.0);
+            double beta = parameters.get<double>("beta",2.0); 
+		    double mu2 = beta*mu1;
+
+	 		auto muTerm = [mu1,mu2,width] (const Domain& z) { 
+		      	if (z[0] < 0 && sqrt(pow(z[1],2) + pow(z[2],2) ) < width/4.0  ||  0 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > width/4.0)
+		        	return mu1;
+		      	else
+		        	return mu2;
+		    };
+            
+            return muTerm;
+   		}
+        else if (imp == "square_fiber"){
+		    double mu1     = parameters.get<double>("mu1",10.0);
+            double beta = parameters.get<double>("beta",2.0); 
+		    double mu2 = beta*mu1;
+
+	 		auto muTerm = [mu1,mu2,width] (const Domain& z) { 
+                if (z[0] < 0 && std::max(abs(z[1]), abs(z[2]))  < width/4.0  ||  0 < z[0] && std::max(abs(z[1]), abs(z[2]))  > width/4.0) 
+		        	return mu1;
+		      	else
+		        	return mu2;
+		    };
+            
+            return muTerm;
+   		}
+        else if (imp == "matrix_material_circles")   // Matrix material with prestrained Fiber inclusions (circular cross-section)
+        {
+            double mu1     = parameters.get<double>("mu1",10.0);
+            double beta    = parameters.get<double>("beta",2.0); 
+		    double mu2     = beta*mu1;
+            
+            int nF    = parameters.get<int>("nF", 2);    //number of Fibers in each Layer
+            double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) );   //default: half of the max-fiber-radius mrF = (width/(2.0*nF)) 
+            
+            
+            
+            assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+            
+            auto muTerm = [mu1,mu2,theta,nF,rF,height,width] (const Domain& x)                
+            {              
+                //TODO if Domain == 1... if Domain == 2 ...
+//                 double domainShift = 1.0/2.0;
+                // TEST
+                double domainShift = 0.0;
+                // shift x to domain [0,1]^3
+                FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+//                 std::cout << "matrixMaterial-MU is used" << std::endl;
+                
+                double output;
+                
+                // determine if point is in upper/lower Layer
+                if ( 0  <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))     // have to evenly space fibers... 
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+//                             
+                            if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+                                output = mu2;
+//                                 return mu2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = mu1;
+//                                 return mu1;
+                        }
+                        else {}
+                    }
+                }
+                else  // lower Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))  
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+                            
+                            if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+                                output = mu2;
+//                                 return mu2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = mu1;
+//                                 return mu1;
+                        }
+                        else{}
+                    }
+
+                }
+                return output;
+            };
+            return muTerm;
+        }
+        else if (imp == "matrix_material_squares")    // Matrix material with prestrained Fiber inclusions (square-shaped cross-section)
+        {
+            double mu1     = parameters.get<double>("mu1",10.0);
+            double beta    = parameters.get<double>("beta",2.0); 
+		    double mu2     = beta*mu1;
+            
+            int nF    = parameters.get<int>("nF", 2);    //number of Fibers in each Layer
+            double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) );   //default: half of the max-fiber-radius mrF = (width/(2.0*nF)) 
+            
+            
+            
+            assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+            
+            auto muTerm = [mu1,mu2,theta,nF,rF,height,width] (const Domain& x)                
+            {              
+                //TODO if Domain == 1... if Domain == 2 ...
+//                 double domainShift = 1.0/2.0;
+                // TEST
+                double domainShift = 0.0;
+                // shift x to domain [0,1]^3
+                FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+//                 std::cout << "matrixMaterial-MU is used" << std::endl;
+                double output;
+                
+                // determine if point is in upper/lower Layer
+                if ( 0  <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))     // have to evenly space fibers... 
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+//                             
+                            if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1])  ) <= rF )
+                                output = mu2;
+//                                 return mu2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = mu1;
+//                                 return mu1;
+                        }
+                    }
+                }
+                else  // lower Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))  
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+                            
+                            if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1])  ) <= rF )
+                                output = mu2;
+//                                 return mu2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = mu1;
+//                                 return mu1;
+                        }
+                    }
+                }
+                return output;
+            };
+            return muTerm;
+        }
+        else if (imp == "bilayer"){	
+		    double mu1     = parameters.get<double>("mu1",10.0);
+            double beta = parameters.get<double>("beta",2.0); 
+		    double mu2 = beta*mu1;
+		    
+		    auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+		      if (isInRotatedPlane(phi, z[dim-2], z[dim-1]))
+		        return mu1;
+		      else
+		        return mu2; 
+            };
+		    
+			return muTerm;
+		}
+		
+
+
+		else if (imp == "helix_poisson"){                                         // interessant!
+		    double mu1     = parameters.get<double>("mu1",10.0);
+            double beta    = parameters.get<double>("beta",2.0); 
+		    double mu2     = beta*mu1;
+
+		    double r  = parameters.get<double>("radius");
+			double r2 = parameters.get<double>("radius_helix");
+
+		    auto muTerm = [mu1, mu2, r, r2] (const Domain& z) {
+		    	std::array<double,2> h = {r*cos(2*M_PI*z[0]), r*sin(2*M_PI*z[0])};
+		      	if ( sqrt(pow(z[1]-h[0],2)+pow(z[2]-h[1],2)) < r2 )                   // distance to the helix?
+		        	return mu1;
+		      	else
+		        	return mu2; 
+		    };
+		    
+			return muTerm;
+		}
+
+		else if (imp == "chess_poisson"){
+		    double mu1     = parameters.get<double>("mu1",10.0);
+            double beta    = parameters.get<double>("beta",2.0); 
+		    double mu2     = beta*mu1;
+
+		    auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+		      	if ( (isInRotatedPlane(phi, 	   z[dim-2], z[dim-1]) and isInRotatedPlane(phi + M_PI/2  , z[dim-2], z[dim-1])) or 
+		      		 (isInRotatedPlane(phi + M_PI, z[dim-2], z[dim-1]) and isInRotatedPlane(phi + 3*M_PI/2, z[dim-2], z[dim-1])) )
+		        	return mu1;
+		      	else
+		        	return mu2; 
+		    };
+		    
+			return muTerm;
+		}
+
+// 		else if (imp == "3Dchess_poisson"){
+// 		    double E1                = parameters.get<double>("E1", 17e6); //Faser
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double mu1 = lameMu(E1, nu1);
+// 
+// 		    double E2                = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double mu2 = lameMu(E2, nu2);
+// 
+// 		    auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// 		    	if (z[0]<0.5)
+// 		      		if ( (isInRotatedPlane(phi, 	   z[1], z[2]) and isInRotatedPlane(phi + M_PI/2  , z[1], z[2])) or 
+// 		      			 (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// 		        		return mu1;
+// 		      		else
+// 		        		return mu2; 
+// 		        else
+// 		        	if ( (isInRotatedPlane(phi,  	   z[1], z[2]) and isInRotatedPlane(phi + M_PI/2  , z[1], z[2])) or 
+// 		      			 (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// 		        		return mu2;
+// 		      		else
+// 		        		return mu1; 
+// 		    };
+// 		    
+// 			return muTerm;
+// 		}
+// 		else if (imp == "vertical_bilayer_poisson"){ 
+//     		double E1                = parameters.get<double>("E1", 17e6); //material1
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double mu1 = lameMu(E1, nu1);
+// 		    double E2                = parameters.get<double>("E2", 17e6); //material2
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double mu2 = lameMu(E2, nu2);
+// 
+// 		    auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// 		      if ( isInRotatedPlane(phi+M_PI/2.0, z[dim-2], z[dim-1]) )//x3
+// 		        return mu1;
+// 		      else
+// 		        return mu2; };
+// 		    
+// 			return muTerm;
+// 		}
+// 
+// 		else if (imp == "microstructured_bilayer_poisson"){ 
+//     		double E1                = parameters.get<double>("E1", 17e6); //material1
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double mu1 = lameMu(E1, nu1);
+// 		    double E2                = parameters.get<double>("E2", 17e6); //material2
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double mu2 = lameMu(E2, nu2);
+// 
+// 		    auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// 		      if ( isInRotatedPlane(phi, z[0]-0.5, z[1]) )                 //TODO understand this
+// 		        return mu1;
+// 		      else
+// 		        return mu2; };
+// 		    
+// 			return muTerm;
+// 		}
+// 
+// 		else { //(imp == "bilayer_poisson")
+//     		double E1                = parameters.get<double>("E1", 17e6); //material1
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double mu1 = lameMu(E1, nu1);
+// 		    double E2                = parameters.get<double>("E2", 17e6); //material2
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double mu2 = lameMu(E2, nu2);
+// 
+// 		    auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// 		      if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) ) //x2
+// 		        return mu1;
+// 		      else
+// 		        return mu2; };
+// 		    
+// 			return muTerm;
+// 		}
+	}
+
+	FuncScalar getLambda(ParameterTree parameters)
+    {
+
+		std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
+		//std::array<std::string,5> imps({"homogen_poisson" "bilayer_poisson" "chess_poisson" "3Dchess_poisson" "vertical_bilayer_poisson"});
+    	//int i_imp =  parameters.get<int>("impnumber");
+    	//std::string imp =  imps[i_imp];
+    	double phi = M_PI*parameters.get<double>("material_angle", 0.0)/180;
+        double theta = parameters.get<double>("theta",1.0/4.0);                        //TODO alle parameter hier laden?
+        
+        double width = parameters.get<double>("width", 1.0);
+        double height = parameters.get<double>("height", 1.0);
+
+		if (imp == "homogenenous"){    
+		    double lambda     = parameters.get<double>("lambda",384.615);
+
+		    auto lambdaTerm = [lambda] (const Domain& z) {return lambda;};
+		    
+			return lambdaTerm;
+		}
+        else if (imp == "isotropic_bilayer")
+    	{
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+    
+ 		   auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+                    if (z[2]>0)                                                   
+                        return lambda1;
+                    else
+                        return lambda2;
+                    };
+           return lambdaTerm;
+		}
+
+        else if (imp == "analytical_Example"){
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+    
+            auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+                
+//                  std::cout << "Analytical-LAMBDA is used" << std::endl;  //TEST
+                    if (abs(z[0]) >= (theta/2.0))
+                        return lambda1;
+                    else
+                        return lambda2;
+                    }; 
+           return lambdaTerm;
+		}
+		else if (imp == "circle_fiber"){
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+
+	 		auto lambdaTerm = [lambda1,lambda2,width] (const Domain& z) // Prestrain inducing twist (bisher nur extension)
+			{ 
+		      	if (z[0] < 0 && sqrt(pow(z[1],2) + pow(z[2],2) ) < width/4.0  ||  0 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > width/4.0)
+		        	return lambda1;
+		      	else
+		        	return lambda2;
+		    };
+            
+           return lambdaTerm;
+   		}
+   		else if (imp == "square_fiber"){
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+
+	 		auto lambdaTerm = [lambda1,lambda2,width] (const Domain& z) // Prestrain inducing twist (bisher nur extension)
+			{ 
+		      	if (z[0] < 0 && std::max(abs(z[1]), abs(z[2]))  < width/4.0  ||  0 < z[0] && std::max(abs(z[1]), abs(z[2]))  > width/4.0) 
+		        	return lambda1;
+		      	else
+		        	return lambda2;
+		    };
+            
+           return lambdaTerm;
+   		}
+        else if (imp == "matrix_material_circles")   // Matrix material with prestrained Fiber inclusions (circular cross-section)
+        {
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+            
+            int nF    = parameters.get<int>("nF", 2);    //number of Fibers in each Layer
+            double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) );   //default: half of the max-fiber-radius mrF = (width/(2.0*nF)) 
+            
+            
+
+            assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+            
+            auto lambdaTerm = [lambda1,lambda2, theta,nF,rF,height,width] (const Domain& x)               
+            {              
+                //TODO if Domain == 1... if Domain == 2 ...
+//                 double domainShift = 1.0/2.0;
+                // TEST
+                double domainShift = 0.0;
+                // shift x to domain [0,1]^3
+                FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+//                 std::cout << "matrixMaterial-MU is used" << std::endl;
+                double output;
+                
+                // determine if point is in upper/lower Layer
+                if ( 0  <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))     // have to evenly space fibers... 
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+//                             
+                            if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+                                output = lambda2;
+//                                 return lambda2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = lambda1;
+//                                 return lambda1;
+                        }
+                    }
+                }
+                else  // lower Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))  
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+                            
+                            if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+                                output = lambda2;
+//                                 return lambda2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = lambda1;
+//                                 return lambda1;
+                        }
+                    }
+
+                }
+                return output;
+            };
+            return lambdaTerm;
+        }
+        else if (imp == "matrix_material_squares")    // Matrix material with prestrained Fiber inclusions (square-shaped cross-section)
+        {
+            double lambda1 = parameters.get<double>("lambda1",0.0);
+            double beta = parameters.get<double>("beta",2.0); 
+            double lambda2 = beta*lambda1;
+            
+            int nF    = parameters.get<int>("nF", 2);    //number of Fibers in each Layer
+            double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) );   //default: half of the max-fiber-radius mrF = (width/(2.0*nF)) 
+            
+            assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+            
+            auto lambdaTerm = [lambda1,lambda2, theta,nF,rF,height,width] (const Domain& x)                
+            {              
+                //TODO if Domain == 1... if Domain == 2 ...
+//                 double domainShift = 1.0/2.0;
+                // TEST
+                double domainShift = 0.0;
+                // shift x to domain [0,1]^3
+                FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+//                 std::cout << "matrixMaterial-MU is used" << std::endl;
+                double output;
+                
+                // determine if point is in upper/lower Layer
+                if ( 0  <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))     // have to evenly space fibers... 
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+//                             
+                            if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1])  ) <= rF )
+                                output = lambda2;
+//                                 return lambda2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = lambda1;
+//                                 return lambda1;
+                        }
+                    }
+                }
+                else  // lower Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/double(nF))*(i+1))  
+                        {
+//                             std::cout << " i : " << i << std::endl;
+//                             printvector(std::cout, s, "s" , "--");
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+//                             printvector(std::cout, Fcenter, "Fcenter" , "--");
+                            
+                            if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1])  ) <= rF )
+                                output = lambda2;
+//                                 return lambda2;                                                         //richtig so? Fiber hat tendenziell größeres mu? 
+                            else
+                                output = lambda1;
+//                                 return lambda1;
+                        }
+                    }
+                }
+                return output;
+            };
+            return lambdaTerm;
+        }
+    	else if (imp == "bilayer_lame"){	
+		    double lambda1     = parameters.get<double>("mu1",384.615);
+		    double lambda2     = parameters.get<double>("mu2",384.615);  
+		    
+		    auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+		      if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) )
+		        return lambda1;
+		      else
+		        return lambda2; };
+		    
+			return lambdaTerm;
+		}
+
+		else if (imp == "helix_poisson"){
+		    double E1                = parameters.get<double>("E1", 17e6); //Faser
+		    double nu1               = parameters.get<double>("nu1", 0.3);
+		    double lambda1 = lameLambda(E1,nu1);
+
+		    double E2                = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+		    double nu2               = parameters.get<double>("nu2", 0.5); 
+		    double lambda2 = lameLambda(E2, nu2);
+
+		    double r = parameters.get<double>("radius");
+			double r2 = parameters.get<double>("radius_helix");
+
+		    auto lambdaTerm = [lambda1, lambda2, r, r2] (const Domain& z) {
+		    	std::array<double,2> h = {r*cos(2*M_PI*z[0]), r*sin(2*M_PI*z[0])};
+		      	MatrixRT ret(0.0);
+		      	if ( sqrt(pow(z[1]-h[0],2)+pow(z[2]-h[1],2)) < r2 )
+		        	return lambda1;
+		      	else
+		        	return lambda2; 
+		    };
+		    
+			return lambdaTerm;
+		}
+
+		else if (imp == "chess_poisson"){
+		    double E1                = parameters.get<double>("E1", 17e6); //Faser
+		    double nu1               = parameters.get<double>("nu1", 0.3);
+		    double lambda1 = lameLambda(E1,nu1);
+
+		    double E2                = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+		    double nu2               = parameters.get<double>("nu2", 0.5); 
+		    double lambda2 = lameLambda(E2, nu2);
+
+		    auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+		      	if ( (isInRotatedPlane(phi, 	   z[dim-2], z[dim-1]) and isInRotatedPlane(phi + M_PI/2  , z[dim-2], z[dim-1])) or 
+		      		 (isInRotatedPlane(phi + M_PI, z[dim-2], z[dim-1]) and isInRotatedPlane(phi + 3*M_PI/2, z[dim-2], z[dim-1])) )
+		        	return lambda1;
+		      	else
+		        	return lambda2; 
+		    };
+		    
+			return lambdaTerm;
+		}
+
+// 		else if (imp == "3Dchess_poisson"){
+// 		    double E1                = parameters.get<double>("E1", 17e6); //Faser
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double lambda1 = lameLambda(E1,nu1);
+// 
+// 		    double E2                = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double lambda2 = lameLambda(E2, nu2);
+// 
+// 		    auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// 		    	if (z[0]<0.5)
+// 		      		if ( (isInRotatedPlane(phi, 	   z[1], z[2]) and isInRotatedPlane(phi + M_PI/2  , z[1], z[2])) or 
+// 		      			 (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// 		        		return lambda1;
+// 		      		else
+// 		        		return lambda2; 
+// 		        else
+// 		        	if ( (isInRotatedPlane(phi,  	   z[1], z[2]) and isInRotatedPlane(phi + M_PI/2  , z[1], z[2])) or 
+// 		      			 (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// 		        		return lambda2;
+// 		      		else
+// 		        		return lambda1; 
+// 		    };
+// 		    
+// 			return lambdaTerm;
+// 		}
+// 
+// 		else if (imp == "vertical_bilayer_poisson"){ 
+//     		double E1                = parameters.get<double>("E1", 17e6); //material1
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double lambda1 = lameLambda(E1,nu1);
+// 		    double E2                = parameters.get<double>("E2", 17e6); //material2
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double lambda2 = lameLambda(E2, nu2);
+// 
+// 		    auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// 		      if ( isInRotatedPlane(phi+M_PI/2.0, z[dim-2], z[dim-1]) )
+// 		        return lambda1;
+// 		      else
+// 		        return lambda2; };
+// 		    
+// 			return lambdaTerm;
+// 		}
+// 
+// 		else if (imp == "microstructured_bilayer_poisson"){ 
+//     		double E1                = parameters.get<double>("E1", 17e6); //material1
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double lambda1 = lameLambda(E1,nu1);
+// 		    double E2                = parameters.get<double>("E2", 17e6); //material2
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double lambda2 = lameLambda(E2, nu2);
+// 
+// 		    auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// 		      if ( isInRotatedPlane(phi, z[0]-0.5, z[1]) )
+// 		        return lambda1;
+// 		      else
+// 		        return lambda2; };
+// 		    
+// 			return lambdaTerm;
+// 		}
+// 
+// 		else { //(imp == "bilayer_poisson")
+//     		double E1                = parameters.get<double>("E1", 17e6); //material1
+// 		    double nu1               = parameters.get<double>("nu1", 0.3);
+// 		    double lambda1 = lameLambda(E1, nu1);
+// 		    double E2                = parameters.get<double>("E2", 17e6); //material2
+// 		    double nu2               = parameters.get<double>("nu2", 0.5); 
+// 		    double lambda2 = lameLambda(E2, nu2);
+// 
+// 		    auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// 		      if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) )
+// 		        return lambda1;
+// 		      else
+// 		        return lambda2; };
+// 		    
+// 			return lambdaTerm;
+// 		}
+	}
+
+
+
+};
+
+
+
+
+
+
+
+
+// TODO  add log here? 
+
+
+template <int dim>
+class PrestrainImp 
+{
+
+public:
+//     static const int dim = 3;
+	static const int dimWorld = 3;
+    using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+//     using GridView = CellGridType::LeafGridView;
+    using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+//     using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+    using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+    using ScalarRT = double;
+    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+    using FuncScalar = std::function< ScalarRT(const Domain&) >;
+    using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+    using FuncVector = std::function< VectorRT(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(std::string imp)
+    Func2Tensor getPrestrain(ParameterTree parameters)
+    {
+        
+        std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
+        
+        double width = parameters.get<double>("width", 1.0);
+        double height = parameters.get<double>("height", 1.0);
+        
+        double theta = parameters.get<double>("theta",1.0/4.0); 
+        double p1 = parameters.get<double>("rho1", 1.0);
+        double alpha = parameters.get<double>("alpha", 2.0);
+        double p2 = alpha*p1;
+        
+//     	if (imp == "isotropic_bilayer")
+//     	{
+//             Func2Tensor B = [p1,p2,theta] (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}};
+//                 else 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: isotropic_bilayer " << std::endl;
+//             return B;
+//         }
+        if (imp == "isotropic_bilayer")
+    	{
+            Func2Tensor B = [p1,p2,theta] (const Domain& x)                 //  ISOTROPIC PRESSURE
+            {             
+                if (x[2] > 0)
+                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                else if (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: isotropic_bilayer " << std::endl;
+            return B;
+        }
+        
+        else if (imp == "analytical_Example")
+        {
+            Func2Tensor B = [p1,theta] (const Domain& x)                // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE   
+            {              
+//                 if (abs(x[0]) < (theta/2) && x[2] < 0 )                             //does not make a difference
+                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: analytical_Example "<< std::endl;
+            return B;
+        }
+        else if (imp == "circle_fiber"){
+
+            Func2Tensor B = [p1,theta,width] (const Domain& x)                // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE   
+            {  
+		      	if (x[0] < 0 && sqrt(pow(x[1],2) + pow(x[2],2) ) < width/4.0  ||  0 < x[0] && sqrt(pow(x[1],2) + pow(x[2],2) ) > width/4.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: circle_fiber"<< std::endl;
+            return B;
+   		}
+        else if (imp == "square_fiber"){
+            
+            Func2Tensor B = [p1,theta,width] (const Domain& x)                // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE   
+            {  
+		      	if (x[0] < 0 && std::max(abs(x[1]), abs(x[2]))  < width/4.0  ||  0 < x[0] && std::max(abs(x[1]), abs(x[2]))  > width/4.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: square_fiber"<< std::endl;
+            return B;
+   		}
+        else if (imp == "matrix_material_circles")   // Matrix material with prestrained Fiber inclusions
+        {
+            int nF    = parameters.get<int>("nF", 3);    //number of Fibers in each Layer
+            double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) );   //default: half of the max-fiber-radius mrF = (width/(2.0*nF)) 
+            
+            assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+            
+            Func2Tensor B = [p1,p2,theta,nF,rF,height,width] (const Domain& x)                
+            {              
+                //TODO if Domain == 1... if Domain == 2 ...
+//                 double domainShift = 1.0/2.0;
+                double domainShift = 0.0;
+                // shift x to domain [0,1]^3
+                FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+                
+                MatrixRT output;
+
+                // determine if point is in upper/lower Layer
+                if ((0.0 <= s[2] && s[2] <= 1.0/2.0)) // upperLayer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/nF)*(i+1))  
+                        {
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+//                             
+                            if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+                                output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+//                                 return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                            else
+                                output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+                        }
+                    }
+                }
+                else  // lower Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/nF)*(i+1))  
+                        {
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0};
+    
+                            if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+                                output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+//                                 return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                            else
+                                output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+                        }
+                    }
+
+                }
+                return output;
+            };
+            std::cout <<" Prestrain Type: matrix_material"<< std::endl;
+            return B;
+        }
+        
+        else if (imp == "matrix_material_squares")   // Matrix material with prestrained Fiber inclusions
+        {
+            int nF    = parameters.get<int>("nF", 3);    //number of Fibers in each Layer
+            double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) );   //default: half of the max-fiber-radius mrF = (width/(2.0*nF)) 
+            
+            assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+            
+            Func2Tensor B = [p1,p2,theta,nF,rF,height,width] (const Domain& x)                
+            {              
+                //TODO if Domain == 1... if Domain == 2 ...
+//                 double domainShift = 1.0/2.0;
+                double domainShift = 0.0;
+                // shift x to domain [0,1]^3
+                FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+                
+                MatrixRT output;
+
+                // determine if point is in upper/lower Layer
+                if ((0.0 <= s[2] && s[2] <= 1.0/2.0)) // upperLayer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/nF)*(i+1))  
+                        {
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+//                             
+                            if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1])  ) <= rF )
+                                output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+//                                 return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                            else
+                                output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+                        }
+                    }
+                }
+                else  // lower Layer
+                {
+                    for(size_t i=0; i<nF ;i++) // running through all the Fibers.. 
+                    {
+                        if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 +  (1.0/nF)*(i+1))  
+                        {
+                            // compute Fiber center
+                            FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0};
+    
+                            if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1])  ) <= rF )
+                                output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+//                                 return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+                            else
+                                output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+//                                 return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+                        }
+                    }
+
+                }
+                return output;
+            };
+            std::cout <<" Prestrain Type: matrix_material"<< std::endl;
+            return B;
+        }
+        else   
+        {
+            // TODO other geometries etc... 
+            
+        }
+        // TODO ANISOTROPIC PRESSURE
+        
+
+    }
+    
+    
+    
+    FuncScalar getNormPrestrain(const FuncMatrix& B)
+    {
+	    FuncScalar normB = [&](const Domain& z) {
+	      	return norm(B(z)); 
+	    };
+	    return normB;  
+    }
+
+
+};
+
+
+
+
+
+
+
+#endif

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

dune/microstructure/vtk_filler.hh

+#ifndef DUNE_MICROSTRUCTURE_VTK_FILLER_HH
+#define DUNE_MICROSTRUCTURE_VTK_FILLER_HH
+
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+//#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/microstructure/prestrain_material_geometry.hh>
+
+
+	using std::string;
+	using namespace Dune;
+	using namespace Functions;
+    using namespace Functions::BasisFactory;
+
+	//static const int dimworld = 3;
+
+	//using GT = UGGrid<dim>; //typename MicrostructuredRodGrid::GT<dim>;
+	//using GV = typename GT::LeafGridView;
+	/*
+	using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+	using VectorRT = FieldVector< double, dimworld>;
+	using ScalarRT = double;
+
+
+	using ScalarCT = std::vector<double>;
+	using VectorCT = BlockVector<FieldVector<double,1> >;
+
+	using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
+*/
+	/*
+	template<class Basis>
+	static auto scalarDiscGlobalBasisFunc(
+		const Basis& feBasis,
+		std::function< double(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& scalarFunc)
+	{
+		ScalarCT values;
+    	Functions::interpolate(feBasis, values, scalarFunc);
+    	auto discGlobalBasisFunc = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(feBasis, values);
+    	return discGlobalBasisFunc;
+	}
+
+	template<class Basis>
+	static auto vectorDiscGlobalBasisFunc(
+		const Basis& feBasis,
+		std::function< FieldVector< double, dimworld>(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& vectorFunc)
+	{
+		VectorCT values;
+    	Functions::interpolate(feBasis, values, vectorFunc);
+    	auto discGlobalBasisFunc = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+    		(Functions::subspaceBasis(feBasis), HierarchicVectorView(values));
+    	return discGlobalBasisFunc;
+	}
+*/
+
+template <int dim>
+class FTKfillerContainer {
+
+public:
+
+	static const int dimworld = 3;
+
+// 	using GT = UGGrid<dim>; //typename MicrostructuredRodGrid::GT<dim>;
+    using GT =YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+	using GV = typename GT::LeafGridView;
+	using Domain = typename GV::template Codim<0>::Geometry::GlobalCoordinate;
+	using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+	using VectorRT = FieldVector< double, dimworld>;
+	using ScalarRT = double;
+	using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+	using FuncVector = std::function< VectorRT(const Domain&) >;
+	using FuncScalar = std::function< ScalarRT(const Domain&) >;
+
+	using ScalarCT = std::vector<double>;
+	using VectorCT = BlockVector<FieldVector<double,1> >;
+	using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
+
+
+    void vtkPrestrainNorm(const GV& gridView, const FuncMatrix& B, string filename)          // Updated
+    {
+	    VTKWriter<typename GT::LeafGridView>  vtkWriter(gridView);
+	    Functions::LagrangeBasis<typename GT::LeafGridView, 1> feBasis(gridView);
+	    FuncScalar normB = PrestrainImp<dim>().getNormPrestrain(B);
+
+// 	    auto discGlobalBasisFunc = scalarDiscGlobalBasisFunc(feBasis, normB);
+        auto discGlobalBasisFunc  = Dune::Functions::makeGridViewFunction(normB, gridView);
+//     auto muLocal = localFunction(muGridF);
+
+
+	    ScalarCT values;
+    	Functions::interpolate(feBasis, values, normB);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(feBasis, values);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+    	vtkWriter.write(filename);
+
+    }
+
+    void vtkPrestrainDirector(const typename GT::LeafGridView& gridView, const FuncVector& b, string filename)
+    {
+    	VTKWriter<typename GT::LeafGridView>  vtkWriter(gridView);
+
+        using namespace Functions::BasisFactory;
+
+    	auto basis_CE = makeBasis(
+	    gridView,
+	    power<dimworld>(
+      	lagrange<1>(),
+      	flatLexicographic()));//flatInterleaved()
+
+    	VectorCT values;
+    	Functions::interpolate(basis_CE, values, b);
+    	auto b_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, values);
+
+      	vtkWriter.addVertexData(b_DGBF, VTK::FieldInfo("prestrain_director", VTK::FieldInfo::Type::vector, dimworld));
+    	vtkWriter.write(filename);
+    }
+
+
+    void vtkMaterialCell(const GV& gridView_CE,
+    	const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkMaterialCell ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_CE);
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addCellData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkMaterialCell done.\n";
+    }
+
+
+ 	  void vtkMaterial3DRod(const GV& gridView_R3D, double eps,
+    	const FuncVector& domain_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblem ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_R3D);
+
+	    // Vector data
+	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+	    VectorCT domain_Coeff;
+	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+
+    void vtkProblemCell(const GV& gridView_CE,
+    	const FuncMatrix& B_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblemCell ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_CE);
+
+	    // Vector data
+	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+    	FuncVector Be1_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[0];
+    	};
+    	VectorCT Be1_Coeff;
+	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+	    FuncVector Be2_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[1];
+    	};
+    	VectorCT Be2_Coeff;
+	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[2];
+    	};
+    	VectorCT Be3_Coeff;
+	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+    void vtkProblemCellBVec(const GV& gridView_CE,
+    	const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblemCellBVec ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_CE);
+
+	    // Vector data
+	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+    	FuncVector Be1_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[0];
+    	};
+    	VectorCT Be1_Coeff;
+	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[1];
+    	};
+    	VectorCT Be2_Coeff;
+	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[2];
+    	};
+    	VectorCT Be3_Coeff;
+	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+    	VectorCT B_vec_Coeff;
+		Functions::interpolate(basis_CE, B_vec_Coeff, B_vec_Func);
+		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, B_vec_Coeff);
+  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+
+    void vtkProblem3DRod(const GV& gridView_R3D, double eps,
+    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblem ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_R3D);
+
+	    // Vector data
+	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+	    VectorCT domain_Coeff;
+	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[0];
+    	};
+    	VectorCT Be1_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[1];
+    	};
+    	VectorCT Be2_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[2];
+    	};
+    	VectorCT Be3_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
+	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+	    {
+	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+	    	return normB_Func(x_Ce);
+	    };
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+     void vtkProblem3DRodBVec(const GV& gridView_R3D, double eps,
+    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblem ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_R3D);
+
+	    // Vector data
+	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+	    VectorCT domain_Coeff;
+	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[0];
+    	};
+    	VectorCT Be1_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[1];
+    	};
+    	VectorCT Be2_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[2];
+    	};
+    	VectorCT Be3_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+	    FuncVector B_vec_R3D_Func = [B_vec_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_vec_Func(x_Ce);
+    	};
+	    VectorCT B_vec_Coeff;
+		Functions::interpolate(basis_R3D, B_vec_Coeff, B_vec_R3D_Func);
+		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_R3D, B_vec_Coeff);
+  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
+	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+	    {
+	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+	    	return normB_Func(x_Ce);
+	    };
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    FuncScalar mu_R3D_Func = [mu_Func, eps](const auto& x)
+	    {
+	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+	    	return mu_Func(x_Ce);
+	    };
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_R3D_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+};
+
+/*
+	template<int dim>
+	void vtkProblemCell(const typename UGGrid<dim>::LeafGridView& gridView_CE,
+    	const std::function< MatrixRT(const typename UGGrid<dim>::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate&)& B_Func,
+    	const std::function< ScalarRT(const typename UGGrid<dim>::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate&)& mu_Func,
+    	std::string filename)
+    {
+    	using GT = UGGrid<dim>;
+    	using GV = typename GT::LeafGridView;
+    	using Domain = typename GV::template Codim<0>::Geometry::GlobalCoordinate;
+    	using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+		using FuncVector = std::function< VectorRT(const Domain&) >;
+		using FuncScalar = std::function< ScalarRT(const Domain&) >;
+    	std::cout << "vtkProblemCell ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_CE);
+
+	    // Vector data
+	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+    	FuncVector Be1_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[0];
+    	};
+    	VectorCT Be1_Coeff;
+	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[1];
+    	};
+    	VectorCT Be2_Coeff;
+	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[2];
+    	};
+    	VectorCT Be3_Coeff;
+	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+	    FuncScalar normB_Func = getNormPrestrain(B_Func);
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+    void vtkProblemCellBVec(const GV& gridView_CE,
+    	const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblemCell ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_CE);
+
+	    // Vector data
+	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+    	FuncVector Be1_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[0];
+    	};
+    	VectorCT Be1_Coeff;
+	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[1];
+    	};
+    	VectorCT Be2_Coeff;
+	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_Func = [B_Func](const auto& x)
+    	{
+    		return B_Func(x)[2];
+    	};
+    	VectorCT Be3_Coeff;
+	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+    	VectorCT B_vec_Coeff;
+		Functions::interpolate(basis_CE, B_vec_Coeff, B_vec_Func);
+		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, B_vec_Coeff);
+  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+	    FuncScalar normB_Func = getNormPrestrain(B_Func);
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+
+    void vtkProblem3DRod(const GV& gridView_R3D, double eps,
+    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblem ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_R3D);
+
+	    // Vector data
+	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+	    VectorCT domain_Coeff;
+	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[0];
+    	};
+    	VectorCT Be1_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[1];
+    	};
+    	VectorCT Be2_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[2];
+    	};
+    	VectorCT Be3_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+	    FuncScalar normB_Func = getNormPrestrain(B_Func);
+	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+	    {
+	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+	    	return normB_Func(x_Ce);
+	    };
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+
+     void vtkProblem3DRodBVec(const GV& gridView_R3D, double eps,
+    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+    	std::string filename)
+    {
+
+    	std::cout << "vtkProblem ...\n";
+
+	    VTKWriter<GV>  vtkWriter(gridView_R3D);
+
+	    // Vector data
+	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+	    VectorCT domain_Coeff;
+	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[0];
+    	};
+    	VectorCT Be1_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[1];
+    	};
+    	VectorCT Be2_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_Func(x_Ce)[2];
+    	};
+    	VectorCT Be3_R3D_Coeff;
+	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+	    FuncVector B_vec_R3D_Func = [B_vec_Func, eps](const auto& x)
+    	{
+    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+    		return B_vec_Func(x_Ce);
+    	};
+	    VectorCT B_vec_Coeff;
+		Functions::interpolate(basis_R3D, B_vec_Coeff, B_vec_R3D_Func);
+		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_R3D, B_vec_Coeff);
+  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+
+	    // Scalar data
+	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+	    FuncScalar normB_Func = getNormPrestrain(B_Func);
+	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+	    {
+	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+	    	return normB_Func(x_Ce);
+	    };
+	    ScalarCT prestrain_Coeff;
+    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, prestrain_Coeff);
+	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+	    FuncScalar mu_R3D_Func = [mu_Func, eps](const auto& x)
+	    {
+	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+	    	return mu_Func(x_Ce);
+	    };
+	    ScalarCT material_Coeff;
+    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_R3D_Func);
+    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+    		(scalarFeBasis, material_Coeff);
+	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+	    vtkWriter.write(filename);
+
+	    std::cout << "vtkProblem done.\n";
+    }
+*/
+
+#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"
+#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.25
+
+
+#############################################
+#  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
+
+
+####  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.25
+#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)
+
+
+
+
+
+
+
+
+)*mu2)
+
+
+
+
+
+
+
+

inputs/computeMuGamma.parset

+# --- Parameter File as Input for 'Compute_MuGamma'
+# 
+# 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"
+#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.25
+
+
+#############################################
+#  Material parameters
+#############################################
+
+write_materialFunctions = true   # VTK mu-functions , lambda-functions
+
+beta=2.0
+
+mu1=10.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
+
+theta=0.2
+
+#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)
+
+
+
+
+
+
+
+
+)*mu2)
+
+
+
+
+
+
+
+

outputs/output.txt

+----- Input Parameters -----: 
+alpha: 5
+gamma: 0.01
+theta: 0.25
+beta: 2
+material parameters: 
+mu1: 10
+mu2: 20
+lambda1: 0
+lambda2: 0
+----------------------------: 
+Number of Elements in each direction: [8,8,8]
+size of FiniteElementBasis: 1728
+Solver-type used:  CG-Solver
+---------- OUTPUT ----------
+ --------------------
+Corrector-Matrix M_1: 
+-4.7875e-09 1.11201e-09 0
+1.11201e-09 -1.65324e-21 0
+0 0 0
+
+ --------------------
+Corrector-Matrix M_2: 
+1.04621e-25 -1.08615e-20 0
+-1.08615e-20 -9.4846e-18 0
+0 0 0
+
+ --------------------
+Corrector-Matrix M_3: 
+1.94872e-09 -9.19575e-10 0
+-9.19575e-10 2.93022e-20 0
+0 0 0
+
+ --------------------
+Effective Matrix Q: 
+2.08099 -5.09371e-24 2.86003e-10
+7.68785e-25 2.08333 -5.64259e-22
+-7.27118e-09 -5.65834e-22 2.08306
+
+--- Prestrain Output --- 
+B_hat: -1.24298 -1.25 2.78718e-08
+B_eff: -0.597302 -0.6 1.12953e-08 (Effective Prestrain)
+------------------------ 
+q1: 2.08099
+q2: 2.08333
+q3: 2.08306
+effective b1: -0.597302
+effective b2: -0.6
+effective b3: 1.12953e-08
+b1_hat: -1.24298
+b2_hat: -1.25
+b3_hat: 2.78718e-08
+mu_gamma=2.08306
+----- print analytic solutions -----
+b1_hat_ana : -0.714286
+b2_hat_ana : -1.25
+b3_hat_ana : 0
+b1_eff_ana : -0.375
+b2_eff_ana : -0.6
+b3_eff_ana : 0
+q1_ana : 1.90476
+q2_ana : 2.08333
+q3 should be between q1 and q2
+   Levels    |  L2SymError  |    Order     |  L2SymNorm   | L2SNorm_ana  |    L2Norm    | 
+------------------------------------------------------------------------------------------
+     3       & 7.09613e-02  & 0.00000e+00  & 8.19784e-03  & 7.14286e-02  & 5.58336e-03  & 
+  Levels     | |q1_ana-q1|  | |q2_ana-q2|  |      q3      | |b1_ana-b1|  | |b2_ana-b2|  | |b3_ana-b3|  | 
+---------------------------------------------------------------------------------------------------------
+     3       & 1.76231e-01  & 8.88178e-15  & 2.08306e+00  & 2.22302e-01  & 2.44249e-15  & 1.12953e-08  & 

src/CMakeLists.txt

+#add_executable("dune-microstructure" dune-microstructure.cc)
+#target_link_dune_default_libraries("dune-microstructure")
+#set(CMAKE_BUILD_TYPE Debug)#
+
+
+set(programs Cell-Problem
+	     Cell-Problem_muGamma
+	     Compute_MuGamma
+             )
+
+
+foreach(_program ${programs})
+  add_executable(${_program} ${_program}.cc)
+  target_link_dune_default_libraries(${_program})
+  #add_dune_pythonlibs_flags(${_program})
+  #target_link_libraries(${_program} tinyxml2)
+#  target_compile_options(${_program} PRIVATE "-fpermissive")
+endforeach()
+
+
+set(CMAKE_BUILD_TYPE Debug)