From e6c12bf00a37cee57935353bcb4fb818f042cbf2 Mon Sep 17 00:00:00 2001 From: romeric Date: Wed, 30 Nov 2016 18:58:58 +0100 Subject: [PATCH] Fixed the factor bug in most electromechanics materials. Two equivalent models one with helmholtz and one with internal energy implemented for proof checks. 3D problems still behave strange, further tests required --- .../BoundaryCondition/BoundaryCondition.py | 2 +- .../LegendreTransform/LegendreTransform.py | 45 ++--- .../IsotropicElectroMechanics_100.py | 32 ++-- .../IsotropicElectroMechanics_101.py | 20 +- .../IsotropicElectroMechanics_102.py | 21 +-- .../IsotropicElectroMechanics_103.py | 21 +-- .../IsotropicElectroMechanics_104.py | 21 +-- .../IsotropicElectroMechanics_105.py | 8 +- .../IsotropicElectroMechanics_106.py | 6 +- .../IsotropicElectroMechanics_107.py | 64 +++---- .../IsotropicElectroMechanics_200.py | 107 +++++++++++ .../IsotropicElectroMechanics_201.py | 104 +++++++++++ .../IsotropicElectroMechanics_3.py | 1 - Florence/MaterialLibrary/Piezoelectric_100.py | 172 ++++++++++-------- .../__LLDispatch__/__init__.py | 0 Florence/MaterialLibrary/__init__.py | 9 +- Florence/Solver/FEMSolver.py | 5 +- Florence/Tensor/Numeric.pyx | 53 +++++- Florence/Tensor/Tensors.py | 132 +++++++------- Florence/Utils/__init__.py | 2 +- Florence/Utils/debug.py | 49 +++-- .../DisplacementPotentialFormulation.py | 4 +- 22 files changed, 566 insertions(+), 312 deletions(-) create mode 100644 Florence/MaterialLibrary/IsotropicElectroMechanics_200.py create mode 100644 Florence/MaterialLibrary/IsotropicElectroMechanics_201.py create mode 100644 Florence/MaterialLibrary/__LLDispatch__/__init__.py diff --git a/Florence/BoundaryCondition/BoundaryCondition.py b/Florence/BoundaryCondition/BoundaryCondition.py index 6c39d2067..7d3266d6a 100644 --- a/Florence/BoundaryCondition/BoundaryCondition.py +++ b/Florence/BoundaryCondition/BoundaryCondition.py @@ -679,7 +679,7 @@ def ApplyDirichletGetReducedMatrices(self, stiffness, F, AppliedDirichlet, LoadF """ # APPLY DIRICHLET BOUNDARY CONDITIONS - for i in range(0,self.columns_out.shape[0]): + for i in range(self.columns_out.shape[0]): F = F - LoadFactor*AppliedDirichlet[i]*stiffness.getcol(self.columns_out[i]) diff --git a/Florence/LegendreTransform/LegendreTransform.py b/Florence/LegendreTransform/LegendreTransform.py index f2e2e276f..869f706dd 100644 --- a/Florence/LegendreTransform/LegendreTransform.py +++ b/Florence/LegendreTransform/LegendreTransform.py @@ -6,7 +6,7 @@ class LegendreTransform(object): def __init__(self, material_internal_energy=None, material_enthalpy=None, - newton_raphson_tolerance=1e-08): + newton_raphson_tolerance=1e-09): self.material_internal_energy = material_internal_energy self.material_enthalpy = material_enthalpy self.newton_raphson_tolerance = newton_raphson_tolerance @@ -19,41 +19,24 @@ def SetIternalEnergy(self, energy): def SetEnthalpy(self, enthalpy): self.material_enthalpy = enthalpy - def InternalEnergyToEnthalpy(self,W_dielectric, W_coupling_tensor, W_elasticity, in_voigt=True): - """Converts the directional derivatives of the free energy of a system to its electric enthalpy. - - in_voigt=True operates in Voigt form inputs - in_voigt=False operates in indicial form inputs - - - Note that the coupling tensor should be symmetric with respect to the first two indices. - Irrespective of the input option (in_voigt), the output is always in Voigt form""" - + def InternalEnergyToEnthalpy(self,W_dielectric, W_coupling_tensor, W_elasticity): + """Converts the directional derivatives of the internal energy of a system to its electric enthalpy. + Note that the transformation needs to be performed first and then the tensors need to be transformed + to Voigt form + """ # DIELECTRIC TENSOR REMAINS THE SAME IN VOIGT AND INDICIAL NOTATION H_dielectric = - np.linalg.inv(W_dielectric) # COMPUTE THE CONSTITUTIVE TENSORS OF ENTHALPY BASED ON THOSE ON INTERNAL ENERGY - if in_voigt: - # IF GIVEN IN VOIGT FORM - H_coupling_tensor = - np.dot(H_dielectric,W_coupling_tensor.T) - H_elasticity = W_elasticity - np.dot(W_coupling_tensor,H_coupling_tensor) - # H_elasticity = W_elasticity + np.dot(W_coupling_tensor,H_coupling_tensor) # not right - H_coupling_tensor = H_coupling_tensor.T - - else: - # IF GIVEN IN INDICIAL FORM - # H_coupling_tensor = - einsum('mi,jkm->jki',H_dielectric,W_coupling_tensor) - # H_elasticity = Voigt(W_elasticity - einsum('ijm,klm->ijkl',W_coupling_tensor,H_coupling_tensor),1) - - H_coupling_tensor = - einsum('ij,kli->klj',H_dielectric,W_coupling_tensor) - H_elasticity = Voigt(W_elasticity - einsum('ijk,klm->ijlm',W_coupling_tensor,einsum('kji',H_coupling_tensor)),1) - - # print(H_coupling_tensor[:,:,0]) - # print(W_coupling_tensor[:,:,0]) - # print("\n") - - H_coupling_tensor = Voigt(H_coupling_tensor,1) + H_coupling_tensor = - einsum('ij,kli->klj',H_dielectric,W_coupling_tensor) + H_elasticity = Voigt(W_elasticity - einsum('ijk,klm->ijlm',W_coupling_tensor,einsum('kji',H_coupling_tensor)),1) + + # print(H_coupling_tensor[:,:,0]) + # print(W_coupling_tensor[:,:,0]) + # print("\n") + + H_coupling_tensor = Voigt(H_coupling_tensor,1) return H_dielectric, H_coupling_tensor, H_elasticity diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_100.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_100.py index 693af23d5..36c7b4f8c 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_100.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_100.py @@ -4,13 +4,12 @@ from .MaterialBase import Material from Florence.LegendreTransform import LegendreTransform from copy import deepcopy -##################################################################################################### - # Simplest electromechanical model with no coupling in terms of internal energy - # W(C,D0) = W_neo(C) + 1/2/eps_1* (D0*D0) -##################################################################################################### class IsotropicElectroMechanics_100(Material): + """Simplest electromechanical model with no coupling in terms of internal energy + W(C,D0) = W_neo(C) + 1/2/eps_1* (D0*D0) + """ def __init__(self, ndim, **kwargs): mtype = type(self).__name__ @@ -20,10 +19,10 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): @@ -40,26 +39,18 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): DD = np.dot(D.T,D)[0,0] self.elasticity_tensor = lamb*(2.*J-1.)*einsum("ij,kl",I,I) + \ - (mu/J - lamb*(J-1))*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + (mu/J - lamb*(J-1))*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) self.coupling_tensor = np.zeros((self.ndim,self.ndim,self.ndim)) self.dielectric_tensor = J/eps_1*np.linalg.inv(b) - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 - # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) @@ -89,7 +80,8 @@ def ElectricDisplacementx(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): # E = ElectricFieldx.reshape(self.ndim,1) # eps_1 = self.eps_1 # D_exact = eps_1/J*np.dot(b,E) - # print np.linalg.norm(D - D_exact) + # # print np.linalg.norm(D - D_exact) # return D_exact + # print D return D \ No newline at end of file diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_101.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_101.py index f8f7e1a42..b392c5620 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_101.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_101.py @@ -19,10 +19,10 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): @@ -46,19 +46,11 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): self.dielectric_tensor = J/eps_1*I - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 - # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_102.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_102.py index c4da1fbba..75877d470 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_102.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_102.py @@ -19,10 +19,10 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): @@ -49,19 +49,12 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): self.dielectric_tensor = 1./eps_1*I - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 - # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) + # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_103.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_103.py index 5c26c4542..7f9318c7d 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_103.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_103.py @@ -19,10 +19,10 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): @@ -46,19 +46,12 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): self.dielectric_tensor = J/eps_1*np.linalg.inv(b) + J/eps_2*I - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 - # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) + # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_104.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_104.py index da7f4e936..e9b7070e6 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_104.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_104.py @@ -19,10 +19,10 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): @@ -51,19 +51,12 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): self.dielectric_tensor = J/eps_1*np.linalg.inv(b) + 1./eps_2*I - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 - # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) + # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_105.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_105.py index dd56ea8db..e270551b4 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_105.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_105.py @@ -51,12 +51,10 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) + # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_106.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_106.py index 760ec2913..ea4057b13 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_106.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_106.py @@ -55,13 +55,11 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): self.dielectric_tensor = J/eps_1*np.linalg.inv(b) + 1./eps_2*I # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # self.coupling_tensor, self.elasticity_tensor, in_voigt=False) E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) + self.coupling_tensor, self.elasticity_tensor) # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_107.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_107.py index 33ac588c4..b3ebcf122 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_107.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_107.py @@ -3,14 +3,14 @@ from Florence.Tensor import trace, Voigt from .MaterialBase import Material from Florence.LegendreTransform import LegendreTransform -##################################################################################################### - # Electromechanical model in terms of internal energy - # W(C,D) = W_mn(C) + 1/2/eps_1 (D0*D0) + 1/2/eps_2/J (FD0*FD0) - # W_mn(C) = u1*C:I+u2*G:I - 2*(u1+2*u2+6*mue)*lnJ + lamb/2*(J-1)**2 + mue*(C:I)**2 -##################################################################################################### class IsotropicElectroMechanics_107(Material): + """Electromechanical model in terms of internal energy + + W(C,D) = W_mn(C) + 1/2/eps_1 (D0*D0) + 1/2/eps_2/J (FD0*FD0) + W_mn(C) = u1*C:I+u2*G:I - 2*(u1+2*u2+6*mue)*lnJ + lamb/2*(J-1)**2 + mue*(C:I)**2 + """ def __init__(self, ndim, **kwargs): mtype = type(self).__name__ @@ -20,10 +20,10 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): @@ -44,6 +44,11 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): DD = np.dot(D.T,D)[0,0] D0D = np.dot(D,D.T) + if self.ndim==2: + trb = trace(b) + 1 + else: + trb = trace(b) + C_mech = 2.*mu2/J*(2*einsum('ij,kl',b,b) - einsum('ik,jl',b,b) - einsum('il,jk',b,b)) +\ 2.*(mu1+2.*mu2+6*mue)/J*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + \ lamb*(2.*J-1.)*einsum("ij,kl",I,I) - lamb*(J-1)*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) +\ @@ -59,26 +64,19 @@ def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): self.elasticity_tensor = C_mech + C_elect + C_reg self.coupling_tensor = 1./eps_2*(einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I) - einsum('ij,k',I,Dx)) + \ - 8*J/eps_e*einsum('ik,j',b,Dx) + 4*J/eps_e*trace(b)*( einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I) ) +\ + 8*J/eps_e*einsum('ik,j',b,Dx) + 4*J/eps_e*trb*( einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I) ) +\ 4.*J**3/mue/eps_e**2*DD* ( einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I) ) + 8.*J**3/mue/eps_e**2*einsum('i,j,k',Dx,Dx,Dx) self.dielectric_tensor = J/eps_1*np.linalg.inv(b) + 1./eps_2*I + \ - 4.*J/eps_e*trace(b)*I + \ + 4.*J/eps_e*trb*I + \ 8.*J**3/mue/eps_e**2*(D0D + 0.5*DD*I) - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 - # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - self.coupling_tensor, self.elasticity_tensor, in_voigt=False) - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) - + self.coupling_tensor, self.elasticity_tensor) + # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) @@ -103,14 +101,19 @@ def CauchyStress(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): DD = np.dot(D.T,D)[0,0] D0D = np.dot(D,D.T) + if self.ndim==2: + trb = trace(b) + 1 + else: + trb = trace(b) + sigma_mech = 2.0*mu1/J*b + \ - 2.0*mu2/J*(trace(b)*b - np.dot(b,b)) -\ + 2.0*mu2/J*(trb*b - np.dot(b,b)) -\ 2.0*(mu1+2*mu2+6*mue)/J*I +\ lamb*(J-1)*I +\ - 4.*mue/J*trace(b)*b + 4.*mue/J*trb*b sigma_electric = 1/eps_2*(D0D - 0.5*DD*I) # REGULARISATION TERMS - sigma_reg = 4.*J/eps_e*(DD*b + trace(b)*D0D) + \ + sigma_reg = 4.*J/eps_e*(DD*b + trb*D0D) + \ 4.0*J**3/mue/eps_e**2*DD*D0D sigma = sigma_mech + sigma_electric + sigma_reg @@ -118,17 +121,4 @@ def CauchyStress(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): def ElectricDisplacementx(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): D = self.legendre_transform.GetElectricDisplacement(self, StrainTensors, ElectricFieldx, elem, gcounter) - - # SANITY CHECK FOR IMPLICIT COMPUTATUTAION OF D - # I = StrainTensors['I'] - # J = StrainTensors['J'][gcounter] - # b = StrainTensors['b'][gcounter] - # E = ElectricFieldx.reshape(self.ndim,1) - # eps_1 = self.eps_1 - # eps_2 = self.eps_2 - # inverse = np.linalg.inv(J/eps_1*np.linalg.inv(b) + 1./eps_2*I) - # D_exact = np.dot(inverse,E) - # print np.linalg.norm(D - D_exact) - # # return D_exact - return D diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_200.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_200.py new file mode 100644 index 000000000..f4ad47f36 --- /dev/null +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_200.py @@ -0,0 +1,107 @@ +import numpy as np +from numpy import einsum +from Florence.Tensor import trace, Voigt +from .MaterialBase import Material +##################################################################################################### + # Simplest Electromechanical Helmoltz energy +##################################################################################################### + + +class IsotropicElectroMechanics_200(Material): + + def __init__(self, ndim, **kwargs): + mtype = type(self).__name__ + super(IsotropicElectroMechanics_200, self).__init__(mtype, ndim, **kwargs) + # REQUIRES SEPARATELY + self.nvar = self.ndim+1 + self.energy_type = "enthalpy" + + # INITIALISE STRAIN TENSORS + from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures + StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") + self.Hessian(StrainTensors,np.zeros((self.ndim,1))) + + def Hessian(self,StrainTensors,ElectricFieldx=0,elem=0,gcounter=0): + + mu1 = self.mu1 + mu2 = self.mu2 + lamb = self.lamb + eps_1 = self.eps_1 + + I = StrainTensors['I'] + J = StrainTensors['J'][gcounter] + b = StrainTensors['b'][gcounter] + + E = ElectricFieldx.reshape(self.ndim,1) + Ex = E.reshape(E.shape[0]) + EE = np.dot(E,E.T) + be = np.dot(b,ElectricFieldx).reshape(self.ndim,1) + + C_mech = 2.*mu2/J*(2*einsum('ij,kl',b,b) - einsum('ik,jl',b,b) - einsum('il,jk',b,b)) +\ + 2.*(mu1+2.*mu2)/J*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + \ + lamb*(2.*J-1.)*einsum("ij,kl",I,I) - lamb*(J-1)*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + + C_elect = -eps_1/J* (einsum("ik,j,l",I,Ex,Ex) + einsum("il,j,k",I,Ex,Ex) + einsum("i,k,jl",Ex,Ex,I) + einsum("i,l,jk",Ex,Ex,I)) + + C_Voigt = Voigt(C_mech + C_elect, 1) + + P_Voigt = eps_1/J*(einsum('ik,j',I,Ex)+einsum('i,jk',Ex,I)) + P_Voigt = Voigt(P_Voigt,1) + + E_Voigt = -eps_1/J*I + + # if elem==0 and gcounter==0: + # from Florence.Tensor import makezero + # makezero(E_Voigt,tol=1e-8) + # print E_Voigt + # exit() + + # Build the Hessian + factor = -1. + H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) + H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) + H_Voigt = np.concatenate((H1,H2),axis=0) + + self.H_VoigtSize = H_Voigt.shape[0] + + return H_Voigt + + + + def CauchyStress(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): + + mu1 = self.mu1 + mu2 = self.mu2 + lamb = self.lamb + eps_1 = self.eps_1 + + I = StrainTensors['I'] + J = StrainTensors['J'][gcounter] + b = StrainTensors['b'][gcounter] + E = ElectricFieldx.reshape(self.ndim,1) + + if self.ndim==2: + trb = trace(b) + 1 + else: + trb = trace(b) + + simga_mech = 2.0*mu1/J*b + \ + 2.0*mu2/J*(trb*b - np.dot(b,b)) -\ + 2.0*(mu1+2*mu2)/J*I +\ + lamb*(J-1)*I + sigma_electric = eps_1/J*np.dot(E,E.T) + sigma = simga_mech + sigma_electric + + + return sigma + + + def ElectricDisplacementx(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): + + J = StrainTensors['J'][gcounter] + b = StrainTensors['b'][gcounter] + E = ElectricFieldx.reshape(self.ndim,1) + eps_1 = self.eps_1 + + D = eps_1/J*E + return D diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_201.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_201.py new file mode 100644 index 000000000..785d238c0 --- /dev/null +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_201.py @@ -0,0 +1,104 @@ +import numpy as np +from numpy import einsum +from Florence.Tensor import trace, Voigt +from .MaterialBase import Material +from Florence.LegendreTransform import LegendreTransform +##################################################################################################### + # Electromechanical model in terms of internal energy + # W(C,D) = W_mn(C) + 1/2/eps_1 (FD0*FD0) + # W_mn(C) = u1*C:I+u2*G:I - 2*(u1+2*u2)*lnJ + lamb/2*(J-1)**2 +##################################################################################################### + + +class IsotropicElectroMechanics_201(Material): + + def __init__(self, ndim, **kwargs): + mtype = type(self).__name__ + super(IsotropicElectroMechanics_201, self).__init__(mtype, ndim, **kwargs) + # REQUIRES SEPARATELY + self.nvar = self.ndim+1 + self.energy_type = "internal_energy" + self.legendre_transform = LegendreTransform() + + if self.ndim == 2: + self.H_VoigtSize = 5 + elif self.ndim == 3: + self.H_VoigtSize = 9 + + def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): + + mu1 = self.mu1 + mu2 = self.mu2 + lamb = self.lamb + eps_1 = self.eps_1 + + I = StrainTensors['I'] + J = StrainTensors['J'][gcounter] + b = StrainTensors['b'][gcounter] + + D = ElectricDisplacementx.reshape(self.ndim,1) + Dx = D.reshape(self.ndim) + DD = np.dot(D.T,D)[0,0] + + self.elasticity_tensor = 2.*mu2/J*(2*einsum('ij,kl',b,b) - einsum('ik,jl',b,b) - einsum('il,jk',b,b)) +\ + 2.*(mu1+2.*mu2)/J*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + \ + lamb*(2.*J-1.)*einsum("ij,kl",I,I) - lamb*(J-1)*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + + self.coupling_tensor = J/eps_1*(einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I)) + + self.dielectric_tensor = J/eps_1*I + + # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART + E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, + self.coupling_tensor, self.elasticity_tensor) + + # BUILD HESSIAN + factor = -1. + H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) + H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) + H_Voigt = np.concatenate((H1,H2),axis=0) + + return H_Voigt + + + + def CauchyStress(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): + + mu1 = self.mu1 + mu2 = self.mu2 + lamb = self.lamb + eps_1 = self.eps_1 + + I = StrainTensors['I'] + J = StrainTensors['J'][gcounter] + b = StrainTensors['b'][gcounter] + D = ElectricDisplacementx.reshape(self.ndim,1) + + if self.ndim==2: + trb = trace(b) + 1 + else: + trb = trace(b) + + simga_mech = 2.0*mu1/J*b + \ + 2.0*mu2/J*(trb*b - np.dot(b,b)) -\ + 2.0*(mu1+2*mu2)/J*I +\ + lamb*(J-1)*I + sigma_electric = J/eps_1*np.dot(D,D.T) + sigma = simga_mech + sigma_electric + + return sigma + + def ElectricDisplacementx(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): + D = self.legendre_transform.GetElectricDisplacement(self, StrainTensors, ElectricFieldx, elem, gcounter) + + + # SANITY CHECK FOR IMPLICIT COMPUTATUTAION OF D + # I = StrainTensors['I'] + # J = StrainTensors['J'][gcounter] + # E = ElectricFieldx.reshape(self.ndim,1) + # eps_1 = self.eps_1 + # D_exact = eps_1/J*E + # # print np.linalg.norm(D - D_exact) + # return D_exact + + return D diff --git a/Florence/MaterialLibrary/IsotropicElectroMechanics_3.py b/Florence/MaterialLibrary/IsotropicElectroMechanics_3.py index 835c40baa..75d5bdbf7 100644 --- a/Florence/MaterialLibrary/IsotropicElectroMechanics_3.py +++ b/Florence/MaterialLibrary/IsotropicElectroMechanics_3.py @@ -22,7 +22,6 @@ def __init__(self, ndim, **kwargs): def Hessian(self,StrainTensors,ElectricFieldx=0,elem=0,gcounter=0): - # Using Einstein summation (using numpy einsum call) d = np.einsum mu = self.mu diff --git a/Florence/MaterialLibrary/Piezoelectric_100.py b/Florence/MaterialLibrary/Piezoelectric_100.py index 1ce9dbc0f..1c59bec5f 100644 --- a/Florence/MaterialLibrary/Piezoelectric_100.py +++ b/Florence/MaterialLibrary/Piezoelectric_100.py @@ -3,6 +3,12 @@ from Florence.Tensor import trace, Voigt from .MaterialBase import Material from Florence.LegendreTransform import LegendreTransform +from math import sqrt +# from numpy import sqrt + +# import pyximport +# pyximport.install(setup_args={'include_dirs': np.get_include()}) +# from Florence.MaterialLibrary.__LLDispatch__._Piezoelectric_100_ import * class Piezoelectric_100(Material): @@ -22,74 +28,75 @@ def __init__(self, ndim, **kwargs): self.energy_type = "internal_energy" self.legendre_transform = LegendreTransform() - self.is_anisotropic = True - - # self.N = np.array([[0,0,1.]]) - self.N = np.array([[0,1.]]).reshape(ndim,1) + self.is_transversely_isotropic = True + if self.ndim==3: + self.H_VoigtSize = 9 + else: + self.H_VoigtSize = 5 - # INITIALISE STRAIN TENSORS - from Florence.FiniteElements.ElementalMatrices.KinematicMeasures import KinematicMeasures - StrainTensors = KinematicMeasures(np.asarray([np.eye(self.ndim,self.ndim)]*2),"nonlinear") - self.Hessian(StrainTensors,np.zeros((self.ndim,1))) def Hessian(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): + # import pyximport + # pyximport.install(setup_args={'include_dirs': np.get_include()}) + # from Florence.MaterialLibrary.__LLDispatch__._Piezoelectric_100_ import _hessian + mu1 = self.mu1 mu2 = self.mu2 mu3 = self.mu3 lamb = self.lamb eps_1 = self.eps_1 eps_2 = self.eps_2 - - N = self.N + eps_3 = self.eps_3 I = StrainTensors['I'] - F = StrainTensors['F'][gcounter] J = StrainTensors['J'][gcounter] b = StrainTensors['b'][gcounter] + F = StrainTensors['F'][gcounter] H = J*np.linalg.inv(F).T - + N = self.anisotropic_orientations[elem][:,None] D = ElectricDisplacementx.reshape(self.ndim,1) + + # H_Voigt = _hessian(F.copy('c'), H.copy('c'), b.copy('c'), J, D, N, mu1, mu2, mu3, lamb, eps_1, eps_2, eps_3) + # self.dielectric_tensor = H_Voigt[3:,3:] + # # print H_Voigt + # return H_Voigt + + FN = np.dot(F,N)[:,0] + HN = np.dot(H,N)[:,0] + innerHN = einsum('i,i',HN,HN) + outerHN = einsum('i,j',HN,HN) Dx = D.reshape(self.ndim) DD = np.dot(D.T,D)[0,0] - HN = np.dot(H,N) - # print N.shape, H.shape, np.dot(H,N)[:,0].shape, np.dot(HN.T,HN).shape - # exit() - HN = np.dot(H,N) - HNHN = np.dot(HN.T,HN)[0,0] - HNOHN = np.dot(HN,HN.T) - HN = HN[:,0] - - C_mech = 2.*mu2/J*(2*einsum('ij,kl',b,b) - einsum('ik,jl',b,b) - einsum('il,jk',b,b)) +\ - 2.*(mu1+2.*mu2)/J*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) + \ - lamb*(2.*J-1.)*einsum("ij,kl",I,I) - lamb*(J-1)*( einsum("ik,jl",I,I)+einsum("il,jk",I,I) ) +\ - 4.0*mu3*(-1./J*einsum('ij,k,l',I,HN,HN) -1./J*einsum('i, j,kl',HN,HN,I) + HNHN * (einsum('ij,kl',I,I) - \ - 0.5*einsum('ik,jl',I,I) - 0.5*einsum('il,jk',I,I) ) + 0.5/J*einsum('il,jk',I,HNOHN) +\ - 0.5/J*einsum('ik,jl',I,HNOHN) + 0.5/J*einsum('jl,ik',I,HNOHN) + 0.5/J*einsum('jk,il',I,HNOHN) ) + # Iso + Aniso + C_mech = 2.*mu2/J* ( 2.0*einsum('ij,kl',b,b) - einsum('ik,jl',b,b) - einsum('il,jk',b,b) ) + \ + 2.*(mu1+2*mu2+mu3)/J * ( einsum('ik,jl',I,I) + einsum('il,jk',I,I) ) + \ + lamb*(2.*J-1.)*einsum('ij,kl',I,I) - lamb*(J-1.) * ( einsum('ik,jl',I,I) + einsum('il,jk',I,I) ) - \ + 4.*mu3/J * ( einsum('ij,kl',I,outerHN) + einsum('ij,kl',outerHN,I) ) + \ + 2.*mu3/J*innerHN*(2.0*einsum('ij,kl',I,I) - einsum('ik,jl',I,I) - einsum('il,jk',I,I) ) + \ + 2.*mu3/J * ( einsum('il,j,k',I,HN,HN) + einsum('jl,i,k',I,HN,HN) + \ + einsum('ik,j,l',I,HN,HN) + einsum('jk,i,l',I,HN,HN) ) C_elect = 1./eps_2*(0.5*DD*(einsum('ik,jl',I,I) + einsum('il,jk',I,I) + einsum('ij,kl',I,I) ) - \ - einsum('ij,k,l',I,Dx,Dx) - einsum('i,j,kl',Dx,Dx,I)) + einsum('ij,k,l',I,Dx,Dx) - einsum('i,j,kl',Dx,Dx,I)) self.elasticity_tensor = C_mech + C_elect - self.coupling_tensor = 1./eps_2*(einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I) - einsum('ij,k',I,Dx)) + + self.coupling_tensor = 1./eps_2*(einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I) - einsum('ij,k',I,Dx)) + \ + 2.*J*sqrt(mu3/eps_3)*(einsum('ik,j',I,Dx) + einsum('i,jk',Dx,I)) + \ + 2.*sqrt(mu3/eps_3)*(einsum('ik,j',I,FN) + einsum('i,jk',FN,I)) - self.dielectric_tensor = J/eps_1*np.linalg.inv(b) + 1./eps_2*I - if self.ndim == 2: - self.H_VoigtSize = 5 - elif self.ndim == 3: - self.H_VoigtSize = 9 + self.dielectric_tensor = J/eps_1*np.linalg.inv(b) + 1./eps_2*I + 2.*J*sqrt(mu3/eps_3)*I # TRANSFORM TENSORS TO THEIR ENTHALPY COUNTERPART - # E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - # self.coupling_tensor, self.elasticity_tensor, in_voigt=False) E_Voigt, P_Voigt, C_Voigt = self.legendre_transform.InternalEnergyToEnthalpy(self.dielectric_tensor, - Voigt(self.coupling_tensor,1), Voigt(self.elasticity_tensor,1), in_voigt=True) + self.coupling_tensor, self.elasticity_tensor) # BUILD HESSIAN - factor = 1. + factor = -1. H1 = np.concatenate((C_Voigt,factor*P_Voigt),axis=1) H2 = np.concatenate((factor*P_Voigt.T,E_Voigt),axis=1) H_Voigt = np.concatenate((H1,H2),axis=0) @@ -104,58 +111,75 @@ def CauchyStress(self,StrainTensors,ElectricDisplacementx,elem=0,gcounter=0): mu2 = self.mu2 mu3 = self.mu3 lamb = self.lamb + eps_1 = self.eps_1 eps_2 = self.eps_2 - - N = self.N + eps_3 = self.eps_3 I = StrainTensors['I'] - F = StrainTensors['F'][gcounter] J = StrainTensors['J'][gcounter] b = StrainTensors['b'][gcounter] + F = StrainTensors['F'][gcounter] H = J*np.linalg.inv(F).T - D = ElectricDisplacementx.reshape(self.ndim,1) + N = self.anisotropic_orientations[elem][:,None] + FN = np.dot(F,N)[:,0] + HN = np.dot(H,N)[:,0] + outerFN = einsum('i,j',FN,FN) + innerHN = einsum('i,i',HN,HN) + outerHN = einsum('i,j',HN,HN) + + D = ElectricDisplacementx.reshape(self.ndim,1) + Dx = D.reshape(self.ndim) DD = np.dot(D.T,D)[0,0] + D0D = np.dot(D,D.T) - FN = np.dot(F,N) - HN = np.dot(H,N) - HNHN = np.dot(HN.T,HN)[0,0] - HNOHN = np.dot(HN,HN.T) - FNOFN = np.dot(FN,FN.T) - HN = HN[:,0] - - simga_mech = 2.0*mu1/J*b + \ - 2.0*mu2/J*(trace(b)*b - np.dot(b,b)) -\ - 2.0*(mu1+2*mu2)/J*I +\ + if self.ndim == 3: + trb = trace(b) + elif self.ndim == 2: + trb = trace(b) + 1 + + sigma_mech = 2.*mu1/J*b + \ + 2.*mu2/J*(trb*b - np.dot(b,b)) - \ + 2.*(mu1+2*mu2+mu3)/J*I + \ lamb*(J-1)*I +\ - 2.0*mu3/J*HNHN*I - 2.*mu3/J*HNOHN + \ - 2.0*mu3/J*FNOFN - sigma_electric = 1/eps_2*(np.dot(D,D.T) - 0.5*DD*I) - sigma = simga_mech + sigma_electric + 2*mu3/J*outerFN +\ + 2*mu3/J*innerHN*I - 2*mu3/J*outerHN - return sigma + sigma_electric = 1./eps_2*(D0D - 0.5*DD*I) +\ + 2.*J*sqrt(mu3/eps_3)*D0D + 2*sqrt(mu3/eps_3)*(einsum('i,j',Dx,FN) + einsum('i,j',FN,Dx)) - def ElectricDisplacementx(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): - D = self.legendre_transform.GetElectricDisplacement(self, StrainTensors, ElectricFieldx, elem, gcounter) + sigma = sigma_mech + sigma_electric - # SANITY CHECK FOR IMPLICIT COMPUTATUTAION OF D - # I = StrainTensors['I'] - # J = StrainTensors['J'][gcounter] - # b = StrainTensors['b'][gcounter] - # E = ElectricFieldx.reshape(self.ndim,1) - # eps_1 = self.eps_1 - # eps_2 = self.eps_2 - # inverse = np.linalg.inv(J/eps_1*np.linalg.inv(b) + 1./eps_2*I) - # D_exact = np.dot(inverse,E) - # print np.linalg.norm(D - D_exact) - # return D_exact + # print sigma + # exit() - return D + return sigma + def ElectricDisplacementx(self,StrainTensors,ElectricFieldx,elem=0,gcounter=0): + # THE ELECTRIC FIELD NEEDS TO BE MODFIED TO TAKE CARE OF CONSTANT TERMS + mu3 = self.mu3 + eps_1 = self.eps_1 + eps_2 = self.eps_2 + eps_3 = self.eps_3 + I = StrainTensors['I'] + J = StrainTensors['J'][gcounter] + b = StrainTensors['b'][gcounter] + F = StrainTensors['F'][gcounter] + H = J*np.linalg.inv(F).T + N = self.anisotropic_orientations[elem][:,None] + FN = np.dot(F,N) + HN = np.dot(H,N) + E = ElectricFieldx.reshape(self.ndim,1) + modElectricFieldx = (E - 2.*sqrt(mu3/eps_3)*FN + 2./J*sqrt(mu3/eps_3)*HN) + + D = self.legendre_transform.GetElectricDisplacement(self, StrainTensors, modElectricFieldx, elem, gcounter) + + # SANITY CHECK FOR IMPLICIT COMPUTATUTAION OF D + # inverse = np.linalg.inv(J/eps_1*np.linalg.inv(b) + 1./eps_2*I + 2.*J*sqrt(mu3/eps_3)*I) + # D_exact = np.dot(inverse, (E - 2.*sqrt(mu3/eps_3)*FN + 2./J*sqrt(mu3/eps_3)*HN) ) + # print np.linalg.norm(D - D_exact) + # return D_exact -# FN*FN = NCN -# S = 2 N O N -# FD0 FD0 = D0 CD0 -# S = 2 D0 O D0 -> sigma = 2J D O D \ No newline at end of file + return D \ No newline at end of file diff --git a/Florence/MaterialLibrary/__LLDispatch__/__init__.py b/Florence/MaterialLibrary/__LLDispatch__/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/Florence/MaterialLibrary/__init__.py b/Florence/MaterialLibrary/__init__.py index 35394ec80..d6956d932 100644 --- a/Florence/MaterialLibrary/__init__.py +++ b/Florence/MaterialLibrary/__init__.py @@ -43,4 +43,11 @@ from IsotropicElectroMechanics_105 import * from IsotropicElectroMechanics_106 import * from IsotropicElectroMechanics_107 import * -from Piezoelectric_100 import * \ No newline at end of file +from Piezoelectric_100 import * + + + + +# HELMHOLTZ AND INTERNAL ENERGY EQUIVALENT MODELS +from IsotropicElectroMechanics_200 import * +from IsotropicElectroMechanics_201 import * \ No newline at end of file diff --git a/Florence/Solver/FEMSolver.py b/Florence/Solver/FEMSolver.py index 64ad3a6d8..44bbf8414 100644 --- a/Florence/Solver/FEMSolver.py +++ b/Florence/Solver/FEMSolver.py @@ -391,7 +391,7 @@ def StaticSolver(self, function_spaces, formulation, solver, K, else: self.norm_residual = np.abs(la.norm(Residual[boundary_condition.columns_in])/self.NormForces) - Eulerx = self.NewtonRaphson(function_spaces, formulation, solver, + Eulerx, Eulerp = self.NewtonRaphson(function_spaces, formulation, solver, Increment,K,NodalForces,Residual,mesh,Eulerx,Eulerp, material,boundary_condition,AppliedDirichletInc) @@ -460,7 +460,6 @@ def NewtonRaphson(self, function_spaces, formulation, solver, Residual[boundary_condition.columns_in] = TractionForces[boundary_condition.columns_in] \ - NodalForces[boundary_condition.columns_in] - # SAVE THE NORM if np.isclose(self.NormForces,0.0): self.norm_residual = np.abs(la.norm(Residual[boundary_condition.columns_in])) @@ -506,7 +505,7 @@ def NewtonRaphson(self, function_spaces, formulation, solver, # break - return Eulerx + return Eulerx, Eulerp diff --git a/Florence/Tensor/Numeric.pyx b/Florence/Tensor/Numeric.pyx index c1abea482..17738637f 100644 --- a/Florence/Tensor/Numeric.pyx +++ b/Florence/Tensor/Numeric.pyx @@ -6,7 +6,7 @@ from libc.math cimport fabs from libcpp.vector cimport vector from cpython cimport bool -__all__ = ['trace','doublecontract','makezero','issymetric','tovoigt', 'cross2d', +__all__ = ['trace','doublecontract','makezero','issymetric','tovoigt', 'tovoigt3', 'cross2d', 'fillin','findfirst','findequal','findequal_approx','findless','findgreater'] # COMPILE USING @@ -209,6 +209,55 @@ cdef _Voigt2(const Real_t *C, Real_t *VoigtA): VoigtA[8] = 0.5*(C[5]+C[6]) + + +@boundscheck(False) +@wraparound(False) +def tovoigt3(np.ndarray[Real_t, ndim=3, mode='c'] e): + """Convert a 3D array to its Voigt represenation""" + cdef np.ndarray[Real_t, ndim=2,mode='c'] VoigtA + cdef int n1dim = e.shape[0] + # DISPATCH CALL TO APPROPRIATE FUNCTION + if n1dim == 3: + VoigtA = np.zeros((6,3),dtype=np.float64) + _Voigt33(&e[0,0,0],&VoigtA[0,0]) + elif n1dim == 2: + VoigtA = np.zeros((3,2),dtype=np.float64) + _Voigt23(&e[0,0,0],&VoigtA[0,0]) + + return VoigtA + + +cdef _Voigt33(const Real_t *e, Real_t *VoigtA): + VoigtA[0] = e[0] + VoigtA[1] = e[9] + VoigtA[2] = e[18] + VoigtA[3] = e[4] + VoigtA[4] = e[13] + VoigtA[5] = e[22] + VoigtA[6] = e[8] + VoigtA[7] = e[17] + VoigtA[8] = e[26] + VoigtA[9] = 0.5*(e[1]+e[3]) + VoigtA[10] = 0.5*(e[10]+e[12]) + VoigtA[11] = 0.5*(e[19]+e[21]) + VoigtA[12] = 0.5*(e[2]+e[6]) + VoigtA[13] = 0.5*(e[11]+e[15]) + VoigtA[14] = 0.5*(e[20]+e[24]) + VoigtA[15] = 0.5*(e[5]+e[7]) + VoigtA[16] = 0.5*(e[14]+e[16]) + VoigtA[17] = 0.5*(e[23]+e[25]) + + +cdef _Voigt23(const Real_t *e, Real_t *VoigtA): + VoigtA[0] = e[0] + VoigtA[1] = e[4] + VoigtA[2] = e[3] + VoigtA[3] = e[7] + VoigtA[4] = 0.5*(e[1]+e[2]) + VoigtA[5] = 0.5*(e[5]+e[6]) + + @boundscheck(False) def cross2d(np.ndarray[double, ndim=2] A, np.ndarray[double, ndim=2] B, str dim="3d", toarray=False): """Cross product of second order tensors A_ij x B_ij defined in the sense of R. de Boer @@ -376,7 +425,7 @@ def findless(np.ndarray A, Real_t num): """An equivalent function to numpy.where(A