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()
""""
Experiment: Taken from
[Rumof,Simon,Smoch]:
TWO-SCALE FINITE ELEMENT APPROXIMATION OF A HOMOGENIZED PLATE MODEL
* Diagonal trusses example
"""
#############################################
# Paths
#############################################
parameterSet.resultPath = '/home/klaus/Desktop/Dune_bendIso/dune-microstructure/outputs_diagonal-trusses'
parameterSet.baseName= 'diagonal-trusses'
#############################################
# Grid parameters
#############################################
nX=8
nY=8
parameterSet.structuredGrid = 'simplex'
parameterSet.lower = '0 0'
parameterSet.upper = '1 1'
parameterSet.elements = str(nX)+' '+ str(nY)
parameterSet.macroGridLevel = 2
#############################################
# 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
parameterSet.vtkwrite_analyticalSolution = 0
parameterSet.vtkWrite_analyticalSurfaceNormal = 0
# 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.VTKwriteMacroPhaseIndicator = True
parameterSet.MacroPhaseSubsamplingRefinement = 3
class GlobalMicrostructure:
""" Class that represents the global microstructure.
The global microstructure can be defined individually
for different (finitely many) macroscopic phases.
Each macroscopic phase corresponds to a 'LocalMicrostructure_'
sub-class that defines the necessary data for the local
micro-problem.
Currently, there are three possibilities for the LocalMicrostructure:
(1) Read the effective quantities (Qhom,Beff) from this parset.
(Constant local microstructure)
(2) Solve the micro-problem once to obtain the effective quantities.
(Constant local microstructure)
(3) Solve for micro-problem for each macroscopic quadrature point.
(Macroscopocally varying local microstructure))
"""
def __init__(self,macroPoint=[0,0]):
self.macroPhases = 1
# define first local microstructure options.
self.macroPhase1_constantMicrostructure = False
self.macroPhase1_readEffectiveQuantities = False;
def macroPhaseIndicator(self,y): #y :macroscopic point
""" Indicatorfunction that determines the domains
i.e. macro-phases of different local microstructures.
"""
return 1;
# Represents the local microstructure in Phase 1
class LocalMicrostructure_1:
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" #hard material
materialRatio = (1.0/50.0)
# self.materialParameters_phase1 = [5.0/3.0, 5.0/2.0]
self.materialParameters_phase1 = [200, 1.0]
#- PHASE 2
self.phase2_type="isotropic"
# self.materialParameters_phase2 = [materialRatio*5.0/3.0, materialRatio*5.0/2.0]
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
#--- Define prestrain function for each phase
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]]
# 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]]
parameterSet.printMicroOutput = 0 # Flag that allows to supress the output of the micro-problem.
#############################################
# Grid parameters
#############################################
parameterSet.microGridLevel = 2
#############################################
# Assembly options
#############################################
parameterSet.set_IntegralZero = 1 #(default = false, This overwrites the option 'set_oneBasisFunction_Zero')
parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
# parameterSet.arbitraryLocalIndex = 7 #(default = 0)
# parameterSet.arbitraryElementNumber = 3 #(default = 0)
parameterSet.cacheElementMatrices = 1 # Use caching of element matrices
#############################################
# Solver Options, Type: #1: CHOLMOD (default) | #2: UMFPACK | #3: GMRES | #4 CG | #5 QR:
# (maybe switch to an iterative solver (e.g. GMRES) if direct solver (CHOLMOD/UMFPACK) runs out of memory for too fine grids)
#############################################
parameterSet.Solvertype = 1
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= 2
# --- Write Correctos to VTK-File:
parameterSet.writeCorrectorsVTK = 0
# --- write effective quantities (Qhom.Beff) to .txt-files
parameterSet.write_EffectiveQuantitiesToTxt = False
#############################################
# Debug options
#############################################
parameterSet.print_debug = 0 #(default=false)
parameterSet.print_conditionNumber = 0
parameterSet.print_corrector_matrices = 0
parameterSet.write_checkOrthogonality = 0 #Check orthogonality (75) from the paper.