add first versions of Experiments by Rumpf et al for validation

experiment/Rumpf-Experiments/create_boundaryBox_twistedvalley.py

+# trace generated using paraview version 5.10.0-RC1
+#import paraview
+#paraview.compatibility.major = 5
+#paraview.compatibility.minor = 10
+
+#### import the simple module from the paraview
+from paraview.simple import *
+#### disable automatic camera reset on 'Show'
+paraview.simple._DisableFirstRenderCameraReset()
+
+# create a new 'Box'
+box1 = Box(registrationName='Box1')
+
+# Properties modified on box1
+box1.XLength = 0.05
+box1.YLength = 1.0
+box1.ZLength = 0.05
+
+# get active view
+renderView1 = GetActiveViewOrCreate('RenderView')
+
+# show data in view
+box1Display = Show(box1, renderView1, 'GeometryRepresentation')
+
+# trace defaults for the display properties.
+box1Display.Representation = 'Surface'
+box1Display.ColorArrayName = [None, '']
+box1Display.SelectTCoordArray = 'TCoords'
+box1Display.SelectNormalArray = 'Normals'
+box1Display.SelectTangentArray = 'None'
+box1Display.OSPRayScaleArray = 'Normals'
+box1Display.OSPRayScaleFunction = 'PiecewiseFunction'
+box1Display.SelectOrientationVectors = 'None'
+box1Display.ScaleFactor = 0.1
+box1Display.SelectScaleArray = 'None'
+box1Display.GlyphType = 'Arrow'
+box1Display.GlyphTableIndexArray = 'None'
+box1Display.GaussianRadius = 0.015700000524520873
+box1Display.SetScaleArray = ['POINTS', 'Normals']
+box1Display.ScaleTransferFunction = 'PiecewiseFunction'
+box1Display.OpacityArray = ['POINTS', 'Normals']
+box1Display.OpacityTransferFunction = 'PiecewiseFunction'
+box1Display.DataAxesGrid = 'GridAxesRepresentation'
+box1Display.PolarAxes = 'PolarAxesRepresentation'
+
+# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
+box1Display.ScaleTransferFunction.Points = [-1.0, 0.0, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0]
+
+# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
+box1Display.OpacityTransferFunction.Points = [-1.0, 0.0, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0]
+
+# reset view to fit data
+renderView1.ResetCamera(False)
+
+# update the view to ensure updated data information
+renderView1.Update()
+
+# Properties modified on box1Display
+box1Display.Position = [0.19, 0.5, 0.0]
+
+# Properties modified on box1Display.DataAxesGrid
+box1Display.DataAxesGrid.Position = [0.19, 0.5, 0.0]
+
+# Properties modified on box1Display.PolarAxes
+box1Display.PolarAxes.Translation = [0.19, 0.5, 0.0]
+
+# change solid color
+# change solid color
+box1Display.AmbientColor = [1.0, 0.8196078431372549, 0.7607843137254902]
+box1Display.DiffuseColor = [1.0, 0.8196078431372549, 0.7607843137254902]
+
+
+
+# create a new 'Box'
+box2 = Box(registrationName='Box2')
+
+# Properties modified on box2
+box2.XLength = 0.05
+box2.YLength = 1.0
+box2.ZLength = 0.05
+
+# get active view
+renderView1 = GetActiveViewOrCreate('RenderView')
+
+# show data in view
+box2Display = Show(box2, renderView1, 'GeometryRepresentation')
+
+# trace defaults for the display properties.
+box2Display.Representation = 'Surface'
+box2Display.ColorArrayName = [None, '']
+box2Display.SelectTCoordArray = 'TCoords'
+box2Display.SelectNormalArray = 'Normals'
+box2Display.SelectTangentArray = 'None'
+box2Display.OSPRayScaleArray = 'Normals'
+box2Display.OSPRayScaleFunction = 'PiecewiseFunction'
+box2Display.SelectOrientationVectors = 'None'
+box2Display.ScaleFactor = 0.3140000104904175
+box2Display.SelectScaleArray = 'None'
+box2Display.GlyphType = 'Arrow'
+box2Display.GlyphTableIndexArray = 'None'
+box2Display.GaussianRadius = 0.015700000524520873
+box2Display.SetScaleArray = ['POINTS', 'Normals']
+box2Display.ScaleTransferFunction = 'PiecewiseFunction'
+box2Display.OpacityArray = ['POINTS', 'Normals']
+box2Display.OpacityTransferFunction = 'PiecewiseFunction'
+box2Display.DataAxesGrid = 'GridAxesRepresentation'
+box2Display.PolarAxes = 'PolarAxesRepresentation'
+
+# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
+box2Display.ScaleTransferFunction.Points = [-1.0, 0.0, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0]
+
+# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
+box2Display.OpacityTransferFunction.Points = [-1.0, 0.0, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0]
+
+# reset view to fit data
+renderView1.ResetCamera(False)
+
+# update the view to ensure updated data information
+renderView1.Update()
+
+# Properties modified on box2Display
+box2Display.Position = [0.81, 0.5, 0.0]
+
+# Properties modified on box2Display.DataAxesGrid
+box2Display.DataAxesGrid.Position = [0.81, 0.5, 0.0]
+
+# Properties modified on box2Display.PolarAxes
+box2Display.PolarAxes.Translation = [0.81, 0.5, 0.0]
+
+# change solid color
+# change solid color
+box2Display.AmbientColor = [1.0, 0.8196078431372549, 0.7607843137254902]
+box2Display.DiffuseColor = [1.0, 0.8196078431372549, 0.7607843137254902]
+
+#================================================================
+# addendum: following script captures some of the application
+# state to faithfully reproduce the visualization during playback
+#================================================================
+
+# get layout
+layout1 = GetLayout()
+
+#--------------------------------
+# saving layout sizes for layouts
+
+# layout/tab size in pixels
+layout1.SetSize(2923, 908)
+
+#-----------------------------------
+# saving camera placements for views
+
+# current camera placement for renderView1
+renderView1.CameraPosition = [0.0, 0.0, 6.07216366868931]
+renderView1.CameraParallelScale = 1.5715916024363863
+
+#--------------------------------------------
+# uncomment the following to render all views
+# RenderAllViews()
+# alternatively, if you want to write images, you can use SaveScreenshot(...).
\ No newline at end of file

experiment/Rumpf-Experiments/diagonal-trusses.py

+import math
+import numpy as np
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+
+
+
+
+#############################################
+#  Paths
+#############################################
+parameterSet.resultPath = '/home/klaus/Desktop/Dune_bendIso/dune-microstructure/outputs_diagonal-trusses'
+parameterSet.outputPath = '/home/klaus/Desktop/Dune_bendIso/dune-microstructure/outputs_diagonal-trusses'
+parameterSet.baseName= 'macrovariationtest'
+
+#############################################
+#  Grid parameters
+#############################################
+nX=8
+nY=8
+
+parameterSet.structuredGrid = 'simplex'
+parameterSet.lower = '0 0'
+parameterSet.upper = '1 1'
+parameterSet.elements = str(nX)+' '+  str(nY)
+
+parameterSet.numLevels = 2   #good results only for refinementLevel >=4
+
+#############################################
+#  Options
+#############################################
+parameterSet.measure_analyticalError = False
+parameterSet.measure_isometryError = False
+parameterSet.vtkWrite_analyticalSurfaceNormal = False
+parameterSet.vtkwrite_analyticalSolution = False
+
+
+parameterSet.conforming_DiscreteJacobian = 0
+#############################################
+#  Solver parameters
+#############################################
+# Tolerance of the multigrid solver
+parameterSet.tolerance = 1e-12
+# Maximum number of multigrid iterations
+parameterSet.maxProximalNewtonSteps = 200
+# Initial regularization
+parameterSet.initialRegularization = 200
+# Measure convergence
+parameterSet.instrumented = 0
+
+############################
+#   Problem specifications
+############################
+# Dimension of the domain (only used for numerical-test python files)
+parameterSet.dim = 2
+
+#############################################
+#  VTK/Output
+#############################################
+# Write discrete solution as .vtk-File
+parameterSet.writeVTK = 1
+
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 0
+
+# Write Dof-Vector to .txt-file
+parameterSet.writeDOFvector = 0
+
+#############################################
+#  Dirichlet boundary indicator
+#############################################
+def dirichlet_indicator(x) :
+    if( (x[0] <= 0.01) or (x[0] >= 0.99)):
+        return True
+    else:
+        return False
+
+
+
+# class IndicatorFunctionClass:
+#     def __init__(self,macroPoint):
+#         self.macroPoint = macroPoint
+
+
+#     def indicator(self,x):
+#         if self.macroPoint >= 5:
+#             print('self =', self.macroPoint)
+#             if (abs(x[0]) < (1.0/4.0) and x[1] <= 0 ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+#         else: 
+#             print('use different indicatorfunction')
+#             if (abs(x[0]) > 2  ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+        
+
+#     def __call__(self,x):           #this defines the operator()
+#         print('using now the operator() ')
+#         if self.macroPoint >= 5:
+#             print('self =', self.macroPoint)
+#             if (abs(x[0]) < (1.0/4.0) and x[1] <= 0 ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+#         else: 
+#             print('use different indicatorfunction')
+#             if (abs(x[0]) > 2  ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+        
+    
+
+
+# class Microstructure:
+#     def __init__(self,macroPoint):
+#         self.macroPoint = macroPoint
+
+
+#     def indicator(self,x):
+#         if self.macroPoint >= 5:
+#             print('self =', self.macroPoint)
+#             if (abs(x[0]) < (1.0/4.0) and x[1] <= 0 ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+#         else: 
+#             print('use different indicatorfunction')
+#             if (abs(x[0]) > 2  ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+        
+
+#     def __call__(self,x):           #this defines the operator()
+#         print('using now the operator() ')
+#         if self.macroPoint >= 5:
+#             print('self =', self.macroPoint)
+#             if (abs(x[0]) < (1.0/4.0) and x[1] <= 0 ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+#         else: 
+#             print('use different indicatorfunction')
+#             if (abs(x[0]) > 2  ):
+#                 return 1    #Phase1   
+#             else :
+#                 return 2    #Phase2
+        
+
+class Microstructure:
+    def __init__(self,macroPoint=[0,0]):  #default value for initialization
+        self.macroPoint = macroPoint
+        # self.macroPoint = macroPoint
+        self.gamma = 1.0    #in the future this might change depending on macroPoint.
+        self.phases = 2     #in the future this might change depending on macroPoint.
+        #--- Define different material phases:
+        #- PHASE 1
+        self.phase1_type="isotropic"
+        self.materialParameters_phase1 = [200, 1.0]   
+        #- PHASE 2
+        self.phase2_type="isotropic"
+        self.materialParameters_phase2 = [100, 1.0]    
+
+
+        self.cacheMacroPhase=True
+        self.macroPhases = [] #store macro phase numbers.
+        self.macroPhaseCount = 0
+
+        self.a = (1.0-macroPoint[0])*(2.0-np.sqrt(3.0))/2.0
+        self.b = macroPoint[0] * (2.0-np.sqrt(3.0))/2.0
+
+        # self.a = (0.5)*(2.0-np.sqrt(3))
+        # self.b = 0.5 * (2.0-np.sqrt(3))
+
+        # self.a = 0
+        # self.b = (2.0-np.sqrt(3))/2
+
+        # print('self.a', self.a)
+        # print('self.b', self.b)
+
+    # --- Indicator function for material phases
+    # def indicatorFunction(self,x):
+    #     if self.macroPoint[0] <= 2.0:
+    #         # if (abs(x[0]) < (theta/2) and x[2] >= 0 ):
+    #         if (abs(x[0]) < (1.0/4.0) and x[2] <= 0 ):
+    #             return 1    #Phase1   
+    #         else :
+    #             return 2    #Phase2
+    #     else:
+    #         if (abs(x[0]) < (1.0/4.0) and x[2] >= 0 ):
+    #             return 1    #Phase1    
+    #         else :
+    #             return 2    #Phase2
+            
+
+    def one_norm(self,v,w):
+        return np.abs(v) + np.abs(w)
+
+
+            
+
+    #--- Indicator function for material phases
+    def indicatorFunction(self,y):
+
+        #shift domain
+        x = [y[0]+0.5,y[1]+0.5]
+        # indicator = (( (self.one_norm(x[0],x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0],x[1])) )\
+        #           | (( (self.one_norm(x[0]-1,x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0]-1,x[1])))\
+        #           | (x[0] < self.a/2.0) | (x[0] >  1-self.a/2.0) | (x[1] < self.a/2.0) | (x[1] >  1-self.a/2.0)  
+    
+        indicator = (( (self.one_norm(x[0],x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0],x[1])) )\
+                | (( (self.one_norm(x[0]-1,x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0]-1,x[1])))\
+                | (x[0] < self.a/2.0) | (x[0] >  1-self.a/2.0) | (x[1] < self.a/2.0) | (x[1] >  1-self.a/2.0)  
+            
+        # indicator = (x[0] < self.b/2.0 )
+        # indicator = (np.abs(x[0]) + np.abs(x[1]) < 1 + self.b/2.0)
+        # indicator = (self.one_norm(x[0],x[1]) < 1 + self.b/2.0)
+
+        if(indicator):
+            return 1    #Phase1   
+        else :
+            return 2    #Phase2
+
+            
+
+
+    def macroPhaseMap(self,y):
+        if y[0] <= 2.0:
+            return 1 #Phase 1
+        else:
+            return 2 #Phase 2
+        
+    
+        
+    #--- Define prestrain function for each phase (also works with non-constant values)
+    def prestrain_phase1(self,x):
+        return [[0, 0, 0], [0,0,0], [0,0,0]]
+        # rho = 0 
+        # # rho = 1.0
+        # return [[rho*1.0, 0, 0], [0,rho*1.0,0], [0,0,rho*1.0]]
+
+    def prestrain_phase2(self,x):
+        return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+
+
+
+
+
+# #Test clamp on other side:
+# def dirichlet_indicator(x) :
+#     if( (x[0] >=3.999) or (x[1]<=0.001)):
+#         return True
+#     else:
+#         return False
+
+
+#boundary-values/derivative function
+# def boundaryValues(x):
+#     # a = 0.0
+#     # a = 0.4
+#     a= 1.4
+#     if(x[0] <= -1.99):
+#         return [x[0] + a, x[1], 0]
+#     if(x[0] >= 1.99):
+#         return [x[0] - a, x[1], 0]
+
+
+# def boundaryValues(x):
+#     return [x[0], x[1], 0]
+
+# def boundaryValuesDerivative(x):
+#     return ((1,0,0),
+#             (0,1,0),
+#             (0,0,1))
+
+#############################################
+#  Microstructure
+############################################
+parameterSet.prestrainFlag = True
+parameterSet.macroscopically_varying_microstructure = False
+parameterSet.read_effectiveQuantities_from_Parset = True # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
+# parameterSet.macroscopically_varying_microstructure = True
+# parameterSet.read_effectiveQuantities_from_Parset = False # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
+
+
+parameterSet.printMicroOutput = False
+
+
+effectivePrestrain= np.array([[0.75, 0.0],
+                              [0.0, 1.0]]);
+
+effectiveQuadraticForm = np.array([[3.0, 0.0, 0.0],
+                                   [0.0, 1.0, 0.0],
+                                   [0.0, 0.0, 0.8]]);
+
+
+# effectivePrestrain= np.array([[1.0, 0.0],
+#                               [0.0, 1.0]]);
+
+# effectiveQuadraticForm = np.array([[1.0, 0.0, 0.0],
+#                                    [0.0, 1.0, 0.0],
+#                                    [0.0, 0.0, 1.0]]);
+#############################################
+#  Initial iterate function
+#############################################
+# def f(x):
+#     return [x[0], x[1], 0]
+
+# def df(x):
+#     return ((1,0),
+#             (0,1),
+#             (0,0))
+# #Rotation:
+# def R(beta):
+#     return  [[math.cos(beta),0],
+#             [0,1],
+#             [-math.sin(beta),0]]
+
+
+# def f(x):
+#     a = 0.5
+#     if(x[0] <= 0.01):
+#         return [x[0], x[1], 0]
+#     elif(x[0] >= 3.99):
+#         return [x[0] - a, x[1], 0]
+#     else:
+#         return [x[0], x[1], 0]
+
+
+# beta = math.pi/4.0
+
+# def df(x):
+#     a = 0.5
+#     if(x[0] <= 0.01):
+#         return R(-1*beta)
+#     elif(x[0] >= 3.99):
+#         return R(beta)
+#     else:
+#         return ((1,0),
+#                 (0,1),
+#                 (0,0))
+    
+
+
+#RUMPF EXAMPLE
+
+
+def f(x):
+    a = (3.0/16.0)
+    if(x[0] <= 0.01):
+        return [x[0]+a, x[1], 0]
+    elif(x[0] >= 0.99):
+        return [x[0] - a, x[1], 0]
+    else:
+        return [x[0], x[1], 0]
+
+
+
+def df(x):
+    return ((1,0),
+            (0,1),
+            (0,0))
+
+
+
+
+
+# def f(x):
+#     # a = 0.4
+#     a = 1.4
+#     if(x[0] <= -1.99):
+#         return [x[0] + a, x[1], 0]
+#     elif(x[0] >= 1.99):
+#         return [x[0] - a, x[1], 0]
+#     else:
+#         return [x[0], x[1], 0]
+
+
+# def df(x):
+#     return ((1,0),
+#             (0,1),
+#             (0,0))
+
+
+
+fdf = (f, df)
+
+
+#############################################
+#  Force
+############################################
+parameterSet.assemble_force_term = True
+
+def force(x):
+    return [0, 0, 1e-2]
+
+
+#############################################
+#  Analytical reference Solution
+#############################################
+# def f(x):
+#     return [x[0], x[1], 0]
+#
+#
+# def df(x):
+#     return ((1,0),
+#             (0,1),
+#             (0,0))
+#
+#
+# fdf = (f, df)
+
+
+
+##################### MICROSCALE PROBLEM ####################
+
+# Microstructure used: Parametrized Laminate.
+
+# --- Choose scale ratio gamma:
+# parameterSet.gamma = 1.0
+# # --- Number of material phases
+# parameterSet.Phases = 4 
+
+# #--- Indicator function for material phases
+# def indicatorFunction(x):
+#     theta=0.25
+#     factor=1
+#     if (abs(x[0]) < (theta/2) and x[2] < 0 ):
+#         return 1    #Phase1
+#     elif (abs(x[0]) > (theta/2) and x[2] > 0 ):
+#         return 2    #Phase2
+#     elif (abs(x[0]) < (theta/2) and x[2] > 0 ):
+#         return 3    #Phase3
+#     else :
+#         return 4    #Phase4
+    
+# ########### Options for material phases: #################################
+# #     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+# #########################################################################
+# ## Notation - Parameter input :
+# # isotropic (Lame parameters) : [mu , lambda]
+# #         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]   # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
+# # transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# # general_anisotropic         : full compliance matrix C
+# ######################################################################
+
+# #--- Define different material phases:
+# #- PHASE 1
+# parameterSet.phase1_type="isotropic"
+# materialParameters_phase1 = [2.0, 0]   
+
+# #- PHASE 2
+# parameterSet.phase2_type="isotropic"
+# materialParameters_phase2 = [1.0, 0]   
+
+# #- PHASE 3
+# parameterSet.phase3_type="isotropic"
+# materialParameters_phase3 = [2.0, 0]
+
+# #- PHASE 4
+# parameterSet.phase4_type="isotropic"
+# materialParameters_phase4 = [1.0, 0]
+
+# #--- Define prestrain function for each phase (also works with non-constant values)
+# def prestrain_phase1(x):
+#     return [[2, 0, 0], [0,2,0], [0,0,2]]
+
+# def prestrain_phase2(x):
+#     return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+# def prestrain_phase3(x):
+#     return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# def prestrain_phase4(x):
+#     return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+
+
+#############################################
+#  Grid parameters
+#############################################
+parameterSet.gridLevel = 2
+
+#############################################
+#  Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1            #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1   #(default = false)
+#parameterSet.arbitraryLocalIndex = 7            #(default = 0)
+#parameterSet.arbitraryElementNumber = 3         #(default = 0)
+
+#############################################
+#  Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 4        # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0  #(default = 2)  degree of information for solver output
+
+
+#############################################
+#  Write/Output options      #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 0   # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1  # VTK norm of B (currently not implemented)
+parameterSet.MaterialSubsamplingRefinement = 1
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0  #(default=false)
+
+parameterSet.print_corrector_matrices = 0
+
+# --- Write Correctos to VTK-File:  
+parameterSet.write_VTK = 0
+
+# --- Use caching of element matrices:  
+parameterSet.cacheElementMatrices = 1
+
+
+
+# --- check orthogonality (75) from paper: 
+parameterSet.write_checkOrthogonality = 0

experiment/Rumpf-Experiments/twisted-valley.py

+import math
+import numpy as np
+
+class ParameterSet(dict):
+    def __init__(self, *args, **kwargs):
+        super(ParameterSet, self).__init__(*args, **kwargs)
+        self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+
+
+
+
+#############################################
+#  Paths
+#############################################
+parameterSet.resultPath = '/home/klaus/Desktop/Dune_bendIso/dune-microstructure/outputs_twisted-valley'
+parameterSet.outputPath = '/home/klaus/Desktop/Dune_bendIso/dune-microstructure/outputs_twisted-valley'
+parameterSet.baseName= 'twisted-valley'
+
+#############################################
+#  Grid parameters
+#############################################
+nX=8
+nY=8
+
+parameterSet.structuredGrid = 'simplex'
+parameterSet.lower = '0 0'
+parameterSet.upper = '1 1'
+parameterSet.elements = str(nX)+' '+  str(nY)
+
+parameterSet.numLevels = 2   #good results only for refinementLevel >=4
+
+#############################################
+#  Options
+#############################################
+parameterSet.measure_analyticalError = False
+parameterSet.measure_isometryError = False
+parameterSet.vtkWrite_analyticalSurfaceNormal = False
+parameterSet.vtkwrite_analyticalSolution = False
+
+
+parameterSet.conforming_DiscreteJacobian = 0
+#############################################
+#  Solver parameters
+#############################################
+# Tolerance of the multigrid solver
+parameterSet.tolerance = 1e-12
+# Maximum number of multigrid iterations
+parameterSet.maxProximalNewtonSteps = 200
+# Initial regularization
+parameterSet.initialRegularization = 200
+# Measure convergence
+parameterSet.instrumented = 0
+
+############################
+#   Problem specifications
+############################
+# Dimension of the domain (only used for numerical-test python files)
+parameterSet.dim = 2
+
+#############################################
+#  VTK/Output
+#############################################
+# Write discrete solution as .vtk-File
+parameterSet.writeVTK = 1
+
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 0
+
+# Write Dof-Vector to .txt-file
+parameterSet.writeDOFvector = 0
+
+#############################################
+#  Dirichlet boundary indicator
+#############################################
+def dirichlet_indicator(x) :
+    if( (x[0] <= 0.01) or (x[0] >= 0.99)):
+        return True
+    else:
+        return False
+
+
+
+
+        
+
+
+#############################################
+#  Microstructure
+############################################
+parameterSet.prestrainFlag = True
+parameterSet.macroscopically_varying_microstructure = False
+parameterSet.read_effectiveQuantities_from_Parset = False # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
+# parameterSet.macroscopically_varying_microstructure = True
+# parameterSet.read_effectiveQuantities_from_Parset = False # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
+
+
+parameterSet.printMicroOutput = True
+
+
+
+
+
+class Microstructure:
+    def __init__(self,macroPoint=[0,0]):  #default value for initialization
+        self.macroPoint = macroPoint
+        # self.macroPoint = macroPoint
+        self.gamma = 0.1    #in the future this might change depending on macroPoint.
+        self.phases = 2     #in the future this might change depending on macroPoint.
+        #--- Define different material phases:
+        #- PHASE 1
+        self.phase1_type="isotropic"
+        self.materialParameters_phase1 = [5.0/3.0, 5.0/2.0]
+
+        self.materialRatio = 1.0/50   
+        #- PHASE 2
+        self.phase2_type="isotropic"
+        self.materialParameters_phase2 = [(5.0/3.0)*self.materialRatio, (5.0/2.0)*self.materialRatio]    
+
+    def one_norm(self,v,w):
+        return np.abs(v) + np.abs(w)
+
+
+    #--- Indicator function for material phases
+    def indicatorFunction(self,y):
+        #shift domain
+        x = [y[0]+0.5,y[1]+0.5]
+
+        indicator = (( (self.one_norm(x[0]-1,x[1])) < (3/4)) & ((self.one_norm(x[0]-1,x[1])) > (1/4)))\
+                | (( (self.one_norm(x[0]-1,x[1])) < (7/4)) & ((self.one_norm(x[0]-1,x[1])) > (5/4)))\
+                
+        # indicator = (( (self.one_norm(x[1]-1,x[0])) < (3/4)) & ((self.one_norm(x[1]-1,x[0])) > (1/4)))\
+        # | (( (self.one_norm(x[1]-1,x[0])) < (7/4)) & ((self.one_norm(x[1]-1,x[0])) > (5/4)))\
+            
+        if(indicator):
+            return 1    #Phase1   
+        else :
+            return 2    #Phase2
+        
+    #--- Define prestrain function for each phase (also works with non-constant values)
+    def prestrain_phase1(self,x):
+        return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+    def prestrain_phase2(self,x):
+        return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+
+
+
+effectivePrestrain= np.array([[0.75, 0.0],
+                              [0.0, 1.0]]);
+
+effectiveQuadraticForm = np.array([[3.0, 0.0, 0.0],
+                                   [0.0, 1.0, 0.0],
+                                   [0.0, 0.0, 0.8]]);
+
+
+# effectivePrestrain= np.array([[1.0, 0.0],
+#                               [0.0, 1.0]]);
+
+# effectiveQuadraticForm = np.array([[1.0, 0.0, 0.0],
+#                                    [0.0, 1.0, 0.0],
+#                                    [0.0, 0.0, 1.0]]);
+#############################################
+#  Initial iterate function
+#############################################
+#RUMPF EXAMPLE
+def f(x):
+    a = (3.0/16.0)
+    if(x[0] <= 0.01):
+        return [x[0]+a, x[1], 0]
+    elif(x[0] >= 0.99):
+        return [x[0] - a, x[1], 0]
+    else:
+        return [x[0], x[1], 0]
+
+
+
+def df(x):
+    return ((1,0),
+            (0,1),
+            (0,0))
+
+
+
+
+
+fdf = (f, df)
+#############################################
+#  Force
+############################################
+parameterSet.assemble_force_term = True
+
+def force(x):
+    return [0, 0, -1e-4]
+
+
+#############################################
+#  Grid parameters
+#############################################
+parameterSet.gridLevel = 4
+
+#############################################
+#  Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1            #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1   #(default = false)
+#parameterSet.arbitraryLocalIndex = 7            #(default = 0)
+#parameterSet.arbitraryElementNumber = 3         #(default = 0)
+
+#############################################
+#  Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 4        # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0  #(default = 2)  degree of information for solver output
+
+
+#############################################
+#  Write/Output options      #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1   # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1  # VTK norm of B (currently not implemented)
+parameterSet.MaterialSubsamplingRefinement = 1
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0  #(default=false)
+
+parameterSet.print_corrector_matrices = 0
+
+# --- Write Correctos to VTK-File:  
+parameterSet.writeCorrectorsVTK = 1
+
+# --- Use caching of element matrices:  
+parameterSet.cacheElementMatrices = 1
+
+
+
+# --- check orthogonality (75) from paper: 
+parameterSet.write_checkOrthogonality = 0