Working coupled phase field with linear models

demos/phase_field_coupled.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Feb  1 08:49:11 2019
+
+@author: bleyerj
+"""
+from dolfin import *
+import mfront_wrapper as mf
+from mfront_wrapper.utils import local_project
+from ufl import Max
+import matplotlib.pyplot as plt
+from mshr import *
+import mgis.behaviour as mgis_bv
+
+N = 50
+mesh = UnitSquareMesh(N, 2*N, "crossed")
+domain = Rectangle(Point(0, 0), Point(1, 1)) - Circle(Point(0.5, 0.5), 0.2, 40)
+mesh = generate_mesh(domain, N)
+plot(mesh, linewidth=0.2)
+plt.show()
+
+Vu = VectorFunctionSpace(mesh, "CG", 1)
+def top(x, on_boundary):
+    return near(x[1], 1) and on_boundary
+def internal(x, on_boundary):
+    return near((x[0]-0.5)**2+(x[1]-0.5)**2, 0.2**2, 0.05) and on_boundary
+Uimp = Expression(("0", "t"), t=0, degree=0)
+bcu = [DirichletBC(Vu, Constant((0, 0)), internal),
+       DirichletBC(Vu, Uimp, top)]
+
+Vd = FunctionSpace(mesh, "CG", 1)
+bcd = DirichletBC(Vd, Constant(0.), internal)
+
+u  = Function(Vu, name="Displacement")
+d = Function(Vd, name="Damage")
+dold = Function(Vd, name="Previous damage")
+
+
+material_u = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
+                                      "TensionCompressionSplit",
+                                      hypothesis="plane_strain",
+                                      material_properties={"YoungModulus": 200.,
+                                                            "PoissonRatio": 0.2})
+problem_u = mf.MFrontNonlinearProblem(u, material_u, quadrature_degree=0)
+problem_u.bc = bcu
+problem_u.integration_type = mgis_bv.IntegrationType.IntegrationWithSecantOperator
+d0 = problem_u.get_state_variable("Damage")
+
+psi =problem_u.get_state_variable("EnergyDensity")
+def solve_u(t):
+    Uimp.t = t
+    d0.interpolate(d)
+    problem_u.affect_state_variables()
+    problem_u.solve(u.vector())
+
+material_d = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
+                                      "PhaseFieldCoupled",
+                                      hypothesis="plane_strain",
+                                      material_properties={"RegularizationLength": 0.02,
+                                                            "FractureEnergy": 1.})
+
+problem_d = mf.MFrontNonlinearProblem(d, material_d, quadrature_degree=0)
+problem_d.register_gradient("Damage", d)
+problem_d.register_gradient("DamageGradient", grad(d))
+H = problem_d.get_state_variable("HistoryFunction")
+# prm = problem_d.solver.parameters
+# prm["report"] = False
+# prm["absolute_tolerance"] = 1e-2
+# prm["relative_tolerance"] = 1e-2
+# prm["maximum_iterations"] = 1
+
+
+
+def solve_d(bc=[]):
+    H.assign(local_project(Max(H, psi), H.function_space(), problem_d.dx))
+    problem_d.affect_state_variables()
+    problem_d.bc = bcd
+    dold.assign(d)
+    problem_d.solve(d.vector())
+    p=plot(d)
+    plt.colorbar(p)
+    plt.show()
+
+tol = 1e-2
+import numpy as np
+for (i, t) in enumerate(np.linspace(0, 2e-1,21)[1:]):
+    res = 1.
+    i = 1
+    print("Time step:", i)
+    while res > tol:
+        solve_d(bcd)
+        solve_u(t)
+        res = np.max(d.vector().get_local()-dold.vector().get_local())
+        print("   Iteration {}: {}".format(i, res))
+        i += 1
+
+

materials/MultiphaseModel.mfront

@@ -8,15 +8,17 @@ e₁.setEntryName("MatrixStrain");
 σ₁.setEntryName("MatrixStress");
 
 @Gradient StrainStensor e₂;
-e₂.setEntryName("InclusionStrain");
+e₂.setEntryName("FiberStrain");
 @Flux StressStensor σ₂;
-σ₂.setEntryName("InclusionStress");
+σ₂.setEntryName("FiberStress");
 
 @Gradient StrainStensor V;
 V.setEntryName("RelativeDisplacement");
 @Flux StressStensor I;
 I.setEntryName("InteractionForce");
 
+//@TangentOperatorBlocks {∂σ₁∕∂Δe₁,∂σ₁∕∂Δe₂,∂σ₂∕∂Δe₁,∂σ₂∕∂Δe₂,∂I∕∂ΔV};
+
 @TangentOperatorBlock ∂σ₁∕∂Δe₁;
 @AdditionalTangentOperatorBlock ∂σ₁∕∂Δe₂;
 @AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₁;
@@ -64,8 +66,8 @@ const Stensor4 Cₕ = {
   ∂σ₁∕∂Δe₂ = λ₁₂ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁₂ ⋅ I₄;
   ∂σ₂∕∂Δe₁ = ∂σ₁∕∂Δe₂;
   ∂σ₂∕∂Δe₂ = C2;
-  σ₁ = ∂σ₁∕∂Δe₁ ⋅ e₁ + ∂σ₁∕∂Δe₂ ⋅ e₂;
-  σ₂ = ∂σ₂∕∂Δe₁ ⋅ e₁ + ∂σ₂∕∂Δe₂ ⋅ e₂;
-  I = kappa*V;
+  σ₁ = ∂σ₁∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₁∕∂Δe₂ ⋅ (e₂ + Δe₂);
+  σ₂ = ∂σ₂∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₂∕∂Δe₂ ⋅ (e₂ + Δe₂);
+  I = kappa*(V + ΔV);
   ∂I∕∂ΔV = kappa ⋅ I₄;
 }

materials/TensionCompressionSplit.mfront

+@Behaviour TensionCompressionSplit;
+@Author Jeremy Bleyer;
+@Date   08/03/2020;
+@Description{
+  "Implementation of phase-field Amor split"
+}
+
+@MaterialProperty stress young;
+young.setGlossaryName("YoungModulus");
+@MaterialProperty real nu;
+nu.setGlossaryName("PoissonRatio");
+
+// Residual stiffness
+@Parameter kres = 1e-4;
+
+@StateVariable real d;
+d.setGlossaryName("Damage");
+
+
+@StateVariable real psi;
+psi.setEntryName("EnergyDensity");
+
+@LocalVariable stress lambda;
+@LocalVariable stress mu;
+@LocalVariable stress kappa;
+@LocalVariable StrainStensor e;
+@LocalVariable StressStensor ed;
+@LocalVariable StressStensor s0;
+@LocalVariable real g;
+
+@InitLocalVariables{
+  lambda = computeLambda(young,nu);
+  mu = computeMu(young,nu);
+  kappa = lambda+2*mu/3.;
+  // degradation function
+  g = pow(1-d, 2)+kres;
+}
+
+@PredictionOperator{
+  if(smt==ELASTIC){
+    Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
+  } else if(smt==SECANTOPERATOR){
+    Dt = g*(lambda*Stensor4::IxI()+2*mu*Stensor4::Id());
+  } else {
+    return FAILURE;
+  }
+}
+
+@Integrator{
+  e  = eto+deto;
+  ed = deviator(e);
+  const strain tr = trace(e);
+  const strain ppe = max(strain(0), tr);
+  const strain pne = min(strain(0), tr);
+  // final stress
+//  sig = g*(kappa*(ppe)*Stensor::Id()+2*mu*ed)+kappa*pne*Stensor::Id();
+  s0 = kappa*tr*Stensor::Id()+2*mu*ed;
+  sig = g*s0;
+  psi = 0.5*(s0|e);
+}
+
+@TangentOperator{
+  if(smt==ELASTIC){
+    Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
+  } else if(smt==SECANTOPERATOR){
+    Dt = g*(lambda*Stensor4::IxI()+2*mu*Stensor4::Id());
+  } else {
+    return FAILURE;
+  }
+}
+

materials/src/Makefile.mfront

@@ -12,7 +12,7 @@ INCLUDES := -I../include \
 
 CXXFLAGS := -Wall -Wfatal-errors -ansi $(shell tfel-config --oflags) -fPIC $(INCLUDES) 
 
-SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiphaseModel-generic.cxx MultiphaseModel.cxx PhaseField-generic.cxx PhaseField.cxx PhaseFieldCoupled-generic.cxx PhaseFieldCoupled.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx
+SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiphaseModel-generic.cxx MultiphaseModel.cxx PhaseField-generic.cxx PhaseField.cxx PhaseFieldCoupled-generic.cxx PhaseFieldCoupled.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx TensionCompressionSplit-generic.cxx TensionCompressionSplit.cxx
 
 makefiles1 = $(SRCCXX:.cxx=.d)
 makefiles2 = $(makefiles1:.cpp=.d)
@@ -22,7 +22,7 @@ makefiles  = $(makefiles2)
 
 all : libBehaviour.so 
 
-libBehaviour.so : MultiphaseModel-generic.o MultiphaseModel.o PhaseField-generic.o PhaseField.o PhaseFieldCoupled-generic.o PhaseFieldCoupled.o IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o StationaryHeatTransfer-generic.o StationaryHeatTransfer.o 
+libBehaviour.so : MultiphaseModel-generic.o MultiphaseModel.o PhaseField-generic.o PhaseField.o PhaseFieldCoupled-generic.o PhaseFieldCoupled.o IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o StationaryHeatTransfer-generic.o StationaryHeatTransfer.o TensionCompressionSplit-generic.o TensionCompressionSplit.o 
 	@$(CXX) -shared $^  -o $@ -L"$(strip $(shell tfel-config --library-path))" $(patsubst %,-l%,$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants)) 
 
 install : 

materials/src/targets.lst

@@ -17,7 +17,9 @@ sources : {
 "LogarithmicStrainPlasticity-generic.cxx",
 "LogarithmicStrainPlasticity.cxx",
 "StationaryHeatTransfer-generic.cxx",
-"StationaryHeatTransfer.cxx"
+"StationaryHeatTransfer.cxx",
+"TensionCompressionSplit-generic.cxx",
+"TensionCompressionSplit.cxx"
 };
 cppflags : {
 "$(shell tfel-config --cppflags --compiler-flags)"
@@ -58,7 +60,12 @@ epts : {
 "StationaryHeatTransfer_Axisymmetrical",
 "StationaryHeatTransfer_PlaneStrain",
 "StationaryHeatTransfer_GeneralisedPlaneStrain",
-"StationaryHeatTransfer_Tridimensional"
+"StationaryHeatTransfer_Tridimensional",
+"TensionCompressionSplit_AxisymmetricalGeneralisedPlaneStrain",
+"TensionCompressionSplit_Axisymmetrical",
+"TensionCompressionSplit_PlaneStrain",
+"TensionCompressionSplit_GeneralisedPlaneStrain",
+"TensionCompressionSplit_Tridimensional"
 };
 };
 headers : {
@@ -85,6 +92,10 @@ headers : {
 "MFront/GenericBehaviour/StationaryHeatTransfer-generic.hxx",
 "TFEL/Material/StationaryHeatTransfer.hxx",
 "TFEL/Material/StationaryHeatTransferBehaviourData.hxx",
-"TFEL/Material/StationaryHeatTransferIntegrationData.hxx"
+"TFEL/Material/StationaryHeatTransferIntegrationData.hxx",
+"MFront/GenericBehaviour/TensionCompressionSplit-generic.hxx",
+"TFEL/Material/TensionCompressionSplit.hxx",
+"TFEL/Material/TensionCompressionSplitBehaviourData.hxx",
+"TFEL/Material/TensionCompressionSplitIntegrationData.hxx"
 };
 };

mfront_wrapper/gradient_flux.py

@@ -51,30 +51,30 @@ class Gradient:
     def update(self, x):
         self.function.vector().set_local(x)
 
-class Grad(Gradient):
-    """ Gradient operator : grad(u) """
-    def __init__(self, variable, name=None):
-        Gradient.__init__(self, variable,
-                          nonsymmetric_tensor_to_vector(grad(variable)), name)
-
-class GradT(Gradient):
-    """ Transposed gradient operator : grad(u).T """
-    def __init__(self, variable, name=None):
-        Gradient.__init__(self, variable,
-                          nonsymmetric_tensor_to_vector(grad(variable).T), name)
-
-class SymGrad(Gradient):
-    """ Symmetric gradient operator : sym(grad(u)) """
-    def __init__(self, variable, name=None):
-        Gradient.__init__(self, variable,
-                          symmetric_tensor_to_vector(sym(grad(variable))), name)
-
-class Grad_Id(Gradient):
-    """ Transformation gradient operator : grad(u) + Id(d) """
-    def __init__(self, variable, name=None):
-        Gradient.__init__(self, variable,
-                          symmetric_tensor_to_vector(grad(variable) +
-                          Identity(variable.geometric_dimension()), T22=1), name)
+# class Grad(Gradient):
+#     """ Gradient operator : grad(u) """
+#     def __init__(self, variable, name=None):
+#         Gradient.__init__(self, variable,
+#                           nonsymmetric_tensor_to_vector(grad(variable)), name)
+
+# class GradT(Gradient):
+#     """ Transposed gradient operator : grad(u).T """
+#     def __init__(self, variable, name=None):
+#         Gradient.__init__(self, variable,
+#                           nonsymmetric_tensor_to_vector(grad(variable).T), name)
+
+# class SymGrad(Gradient):
+#     """ Symmetric gradient operator : sym(grad(u)) """
+#     def __init__(self, variable, name=None):
+#         Gradient.__init__(self, variable,
+#                           symmetric_tensor_to_vector(sym(grad(variable))), name)
+
+# class Grad_Id(Gradient):
+#     """ Transformation gradient operator : grad(u) + Id(d) """
+#     def __init__(self, variable, name=None):
+#         Gradient.__init__(self, variable,
+#                           symmetric_tensor_to_vector(grad(variable) +
+#                           Identity(variable.geometric_dimension()), T22=1), name)
 class Var(Gradient):
     """ The variable itself : u """
     def __init__(self, variable, name=None):

mfront_wrapper/nonlinear_problem.py

@@ -193,12 +193,6 @@ class MFrontNonlinearProblem(NonlinearProblem):
         self.block_shapes = self.material.get_tangent_block_sizes()
         self.flattened_block_shapes = [s[0]*s[1] if len(s)==2 else s[0] for s in self.block_shapes]
         self.update_tangent_blocks()
-        buff = 0
-        for (i, s) in enumerate(self.material.get_state_variable_names()):
-            state_var = self.state_variables[s]
-            block_shape = self.material.get_state_variable_sizes()[i]
-            self.material.data_manager.s0.internal_state_variables[:,buff:buff+block_shape] = state_var.vector().get_local().reshape(self.material.data_manager.n, block_shape)
-            buff += block_shape
 
     def update_tangent_blocks(self):
         buff = 0
@@ -257,6 +251,14 @@ class MFrontNonlinearProblem(NonlinearProblem):
             self.material.data_manager.s0.gradients[:, buff:buff+block_shape] = grad_vals
             buff += block_shape
 
+    def affect_state_variables(self):
+        buff = 0
+        for (i, s) in enumerate(self.material.get_state_variable_names()):
+            state_var = self.state_variables[s]
+            block_shape = self.material.get_state_variable_sizes()[i]
+            self.material.data_manager.s0.internal_state_variables[:,buff:buff+block_shape] = state_var.vector().get_local().reshape(self.material.data_manager.n, block_shape)
+            buff += block_shape
+
     def update_state_variables(self):
         buff = 0
         for (i, s) in enumerate(self.material.get_state_variable_names()):