import math
from python_matrix_operations import *
import ctypes
import os
import sys


Phases=3

mu_ = [80, 80, 60]
lambda_ = [80, 80, 25]


phase1_type="isotropic"
materialParameters_phase1 = [80, 80]

phase2_type="isotropic"
materialParameters_phase2 = [80, 80]

phase3_type="isotropic"
materialParameters_phase3 = [60, 25]


# print('mu_:', mu_)
# A = [[1, 5, 0], [5,1,0], [5,0,1]]
#
# print("sym(A):", sym(A))


# dir_path = os.path.dirname(os.path.realpath("/home/klaus/Desktop/Dune-Testing/dune-microstructure/dune/microstructure/matrix_operations.hh"))
# handle = ctypes.CDLL(dir_path)
#
# handle.create2Darray.argtypes = [ctypes.c_int, ctypes.c_double, ctypes.c_double]
#
# def create2Darray(nside, mx, my):
#     return handle.create2Darray(nside, mx, my)


#Indicator function that determines both phases
# x[0] : y1-component
# x[1] : y2-component
# x[2] : x3-component
#    --- replace with your definition of indicatorFunction:
###############
# Wood
###############
# def f(x):
#     theta=0.25
#     # mu_ = [3, 5, 6]
#     # lambda_ = [9, 7, 8]
#     # mu_ = 3 5 6
#     # lambda_ = 9 7 8

#     if ((abs(x[0]) < theta/2) and x[2]<0.25):
#         return [mu_[0], lambda_[0]]    #latewood
#         # return 5    #latewood
#     elif ((abs(x[0]) > theta/2) and x[2]<0.25):
#         return [mu_[1], lambda_[1]]    #latewood
#         # return 2
#     else :
#         return [mu_[2],lambda_[2]]    #latewood    #Phase3
#         # return 1

def indicatorFunction(x):
    theta=0.25
    factor=1
    if (x[0] <-1/2+theta and x[2]<-1/2+theta):
        return 1    #Phase1
    elif (x[1]< -1/2+theta and x[2]>1/2-theta):
        return 2    #Phase2
    else :
        return 0    #Phase3


def L1(G):
    return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id()   #Phase1

def L2(G):
    return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id()   #Phase1

def L3(G):
    return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id()   #Phase1


# TEST

# def L1(G):
#     return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) ))  #Phase1

# def L2(G):
#     return Add(smult(2.0 * mu_[1], sym(G)),smult(lambda_[1] ,smult(trace(sym(G)),Id()) ))   #Phase1

# def L3(G):
#     return Add(smult(2.0 * mu_[2], sym(G)),smult(lambda_[2] ,smult(trace(sym(G)),Id()) ))   #Phase1


#Workaround
def muValue(x):
    return mu_

def lambdaValue(x):
    return lambda_





###############
# Cross
###############
# def f(x):
#     theta=0.25
#     factor=1
#     if (x[0] <-1/2+theta and x[2]<-1/2+theta):
#         return 1    #Phase1
#     elif (x[1]< -1/2+theta and x[2]>1/2-theta):
#         return 2    #Phase2
#     else :
#         return 0    #Phase3





# def f(x):
#     # --- replace with your definition of indicatorFunction:
#     if (abs(x[0]) > 0.25):
#         return 1    #Phase1
#     else :
#         return 0    #Phase2

def b1(x):
    return [[1, 0, 0], [0,1,0], [0,0,1]]

def b2(x):
    return [[1, 0, 0], [0,1,0], [0,0,1]]

def b3(x):
    return [[0, 0, 0], [0,0,0], [0,0,0]]

# mu=80 80 60

# lambda=80 80 25


# --- elasticity tensor
# def L(G,x):
# def L(G):
#     # input:  symmetric matrix G, position x
#     # output: symmetric matrix LG
#     return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]









# --- elasticity tensor
def L(G,x):
    # input:  symmetric matrix G, position x
    # output: symmetric matrix LG
    theta=0.25
    if (x[0] <-1/2+theta and x[2]<-1/2+theta):
        return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id()   #Phase1
    elif (x[1]< -1/2+theta and x[2]>1/2-theta):
        return 2.0 * mu_[1] * sym(G) + lambda_[1] * trace(sym(G)) * Id()   #Phase2
    else :
        return 2.0 * mu_[2] * sym(G) + lambda_[2] * trace(sym(G)) * Id()   #Phase3
# #     # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();


# def L(G,x):
#     # input:  symmetric matrix G, position x
#     # output: symmetric matrix LG
#     theta=0.25
#     if (x[0] <-1/2+theta and x[2]<-1/2+theta):
#         # return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id()   #Phase1
#         return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) ))
#     elif (x[1]< -1/2+theta and x[2]>1/2-theta):
#         return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) ))   #Phase2
#     else :
#         return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) ))   #Phase3
#         # return [[0, 0, 0], [0,0,0], [0,0,0]]  #Phase3


##TEST
# def L(G,x):
#     # input:  symmetric matrix G, position x
#     # output: symmetric matrix LG
#     theta=0.25
#     if (x[0] <-1/2+theta and x[2]<-1/2+theta):
#         return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
#     elif (x[1]< -1/2+theta and x[2]>1/2-theta):
#         return [[x[0], 1, x[0]], [1, 1, 1],[x[0], x[0], x[0]]]
#     else :
#         return [[0, x[2], x[2]], [0,x[2],0], [0,0,0]]


##TEST
# def L(G,x):
#     # input:  symmetric matrix G, position x
#     # output: symmetric matrix LG
#     theta=0.25
#     if (x[0] <-1/2+theta and x[2]<-1/2+theta):
#         return sym([[1, 1, 1], [1, 1, 1],[1, 1, 1]])
#     elif (x[1]< -1/2+theta and x[2]>1/2-theta):
#         return sym([[x[0], 1, x[0]], [1, 1, 1],[x[0], x[0], x[0]]])
#     else :
#         return sym([[0, x[2], x[2]], [0,x[2],0], [0,0,0]])



# # small speedup..
# def L(G,x):
#     # input:  symmetric matrix G, position x
#     # output: symmetric matrix LG
#     theta=0.25
#     if (x[0] <-1/2+theta and x[2]<-1/2+theta):
#         return mu_[0] * (np.array(G).transpose() + np.array(G)) + lambda_[0] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3)   #Phase1
#     elif (x[1]< -1/2+theta and x[2]>1/2-theta):
#         return mu_[1] * (np.array(G).transpose() + np.array(G)) + lambda_[1] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3)   #Phase2
#     else :
#         return mu_[2] * (np.array(G).transpose() + np.array(G)) + lambda_[2] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3)   #Phase3
#     # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();






# def H(G,x):
#     # input:  symmetric matrix G, position x
#     # output: symmetric matrix LG
#     if (abs(x[0]) > 0.25):
#         return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
#     else:
#         return [[0, 0, 0], [0,0,0], [0,0,0]]


def H(G,x):
    # input:  symmetric matrix G, position x
    # output: symmetric matrix LG
    if (abs(x[0]) > 0.25):
        return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
    else:
        return [[0, 0, 0], [0,0,0], [0,0,0]]

# 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();