Fixed the factor bug in most electromechanics materials. Two equivale…

…nt models one with helmholtz and one with internal energy implemented for proof checks. 3D problems still behave strange, further tests required

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])
 
 

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

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

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)

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)

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)

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)

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)

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)