Adds phase field demo

demos/multiphase_model.py

@@ -64,18 +64,5 @@ plt.figure()
 plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "ob", label=r"$u_x$ (matrix, size={})".format(mat_prop["Size"]))
 plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "or", label=r"$u_x$ (fiber, size={})".format(mat_prop["Size"]))
 
-# mat_prop["Size"] = 0.02
-
-# problem2 = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
-# problem2.bc = bc
-# problem2.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True)
-# problem2.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True)
-# problem2.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True)
-# u = Function(V)
-# problem2.solve(u.vector())
-
-# plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "sb", markerfacecolor="None", label=r"$u_x$ (matrix, size={})".format(mat_prop["Size"]))
-# plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "sr", markerfacecolor="None", label=r"$u_x$ (fiber, size={})".format(mat_prop["Size"]))
-
 plt.legend()
 plt.show()
\ No newline at end of file

demos/phase_field.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
+
+N = 100
+mesh = UnitSquareMesh(N, 2*N, "crossed")
+
+
+V = FunctionSpace(mesh, "CG", 1)
+def crack(x, on_boundary):
+    return near(x[0], 0.5, 0.2) and near(x[1], 0.5)
+def border(x, on_boundary):
+    return on_boundary
+
+bc = DirichletBC(V, Constant(1.), crack)
+
+d = Function(V, name="Damage")
+
+material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
+                                      "PhaseField",
+                                      hypothesis="plane_strain",
+                                      material_properties={"RegularizationLength": 0.02})
+
+problem = mf.MFrontNonlinearProblem(d, material, quadrature_degree=0)
+problem.bc = bc
+problem.register_gradient("Damage", d)
+problem.register_gradient("DamageGradient", grad(d))
+problem.solve(d.vector())
+
+import matplotlib.pyplot as plt
+p=plot(d)
+plt.colorbar(p)
+plt.show()

materials/MultiphaseModel.mfront

 @DSL DefaultGenericBehaviour;
 @Behaviour MultiphaseModel;
-@ModellingHypotheses {PlaneStrain};
+@ModellingHypotheses {PlaneStrain, Tridimensional};
 
 @Gradient StrainStensor e₁;
 e₁.setEntryName("MatrixStrain");
@@ -40,6 +40,14 @@ s.setEntryName("Size");
 @Integrator {
   // remove useless warnings, as we always compute the tangent operator
   static_cast<void>(computeTangentOperator_);
+const Stensor4 Cₕ = {
+      1, 2, 3, 4, 5, 6, // 1ère ligne
+      5, 6, 7, 8, 9, 10, // 2ème ligne
+      1, 2, 3, 4, 11, 12, // 3ème ligne
+      1, 2, 3, 4, 13, 14, // 4ème ligne
+      1, 2, 3, 4, 13, 14, // 4ème ligne
+      1, 2, 3, 4, 13, 14 // 4ème ligne
+  };
   const auto λ₁ = computeLambda(Y1,nu1);
   const auto μ₁ = computeMu(Y1,nu1);
   const auto λ₂ = rho*computeLambda(Y2,nu2);
@@ -47,10 +55,15 @@ s.setEntryName("Size");
   constexpr const auto λ₁₂ = 0.;
   constexpr const auto μ₁₂ = 0.;
   const auto kappa = μ₁/12/(1-rho)/pow(s,2);
-  ∂σ₁∕∂Δe₁ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄;
+  const auto C1 = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄;
+  const auto C2 = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄;
+  const auto idC = invert(C2-C1);
+  const auto avgC = rho*C1+(1-rho)*C2;
+  const auto D11 = (1-rho)*C1-(C1*idC*(Cₕ-avgC)*idC*C1);
+  ∂σ₁∕∂Δe₁ = C1;
   ∂σ₁∕∂Δe₂ = λ₁₂ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁₂ ⋅ I₄;
   ∂σ₂∕∂Δe₁ = ∂σ₁∕∂Δe₂;
-  ∂σ₂∕∂Δe₂ = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄;
+  ∂σ₂∕∂Δe₂ = C2;
   σ₁ = ∂σ₁∕∂Δe₁ ⋅ e₁ + ∂σ₁∕∂Δe₂ ⋅ e₂;
   σ₂ = ∂σ₂∕∂Δe₁ ⋅ e₁ + ∂σ₂∕∂Δe₂ ⋅ e₂;
   I = kappa*V;

materials/PhaseField.mfront

+@DSL DefaultGenericBehaviour;
+@Behaviour PhaseField;
+
+@Gradient TVector g;
+g.setEntryName("DamageGradient");
+@Flux TVector q;
+q.setEntryName("DualDamageGradient");
+
+@Gradient real d;
+d.setGlossaryName("Damage");
+@Flux real Y;
+Y.setEntryName("EnergyRelease");
+
+@TangentOperatorBlock ∂q∕∂Δg;
+@AdditionalTangentOperatorBlock ∂Y∕∂Δd;
+
+@MaterialProperty real l0;
+l0.setEntryName("RegularizationLength");
+
+@ProvidesTangentOperator;
+@Integrator {
+  // remove useless warnings, as we always compute the tangent operator
+  static_cast<void>(computeTangentOperator_);
+  ∂q∕∂Δg = pow(l0, 2)* tmatrix<N, N, real>::Id();
+  ∂Y∕∂Δd = 1;
+  q = pow(l0, 2)*g;
+  Y = ∂Y∕∂Δd ⋅ d;
+}

materials/src/Makefile.mfront

@@ -12,7 +12,7 @@ INCLUDES := -I../include \
 
 CXXFLAGS := -Wall -Wfatal-errors -ansi $(shell tfel-config --oflags) -fPIC $(INCLUDES) 
 
-SRCCXX = MultiphaseModel-generic.cxx MultiphaseModel.cxx
+SRCCXX = MultiphaseModel-generic.cxx MultiphaseModel.cxx PhaseField-generic.cxx PhaseField.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 
+libBehaviour.so : MultiphaseModel-generic.o MultiphaseModel.o PhaseField-generic.o PhaseField.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

@@ -7,7 +7,9 @@ suffix : "so";
 install_path : "";
 sources : {
 "MultiphaseModel-generic.cxx",
-"MultiphaseModel.cxx"
+"MultiphaseModel.cxx",
+"PhaseField-generic.cxx",
+"PhaseField.cxx"
 };
 cppflags : {
 "$(shell tfel-config --cppflags --compiler-flags)"
@@ -22,13 +24,23 @@ link_libraries : {
 "$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants)"
 };
 epts : {
-"MultiphaseModel_PlaneStrain"
+"MultiphaseModel_PlaneStrain",
+"MultiphaseModel_Tridimensional",
+"PhaseField_AxisymmetricalGeneralisedPlaneStrain",
+"PhaseField_Axisymmetrical",
+"PhaseField_PlaneStrain",
+"PhaseField_GeneralisedPlaneStrain",
+"PhaseField_Tridimensional"
 };
 };
 headers : {
 "MFront/GenericBehaviour/MultiphaseModel-generic.hxx",
 "TFEL/Material/MultiphaseModel.hxx",
 "TFEL/Material/MultiphaseModelBehaviourData.hxx",
-"TFEL/Material/MultiphaseModelIntegrationData.hxx"
+"TFEL/Material/MultiphaseModelIntegrationData.hxx",
+"MFront/GenericBehaviour/PhaseField-generic.hxx",
+"TFEL/Material/PhaseField.hxx",
+"TFEL/Material/PhaseFieldBehaviourData.hxx",
+"TFEL/Material/PhaseFieldIntegrationData.hxx"
 };
 };

mfront_wrapper/gradient_flux.py

@@ -23,7 +23,7 @@ class Gradient:
         if len(shape)==1:
             self.shape = shape[0]
         elif shape==():
-            self.shape = 0
+            self.shape = 1
             #self.expression = as_vector([self.expression])
         else:
             self.shape = shape

mfront_wrapper/nonlinear_problem.py

@@ -2,6 +2,7 @@ from dolfin import *
 from mfront_wrapper.utils import *
 from mfront_wrapper.gradient_flux import *
 import mgis.behaviour as mgis_bv
+from ufl import shape
 
 class MFrontNonlinearProblem(NonlinearProblem):
     def __init__(self, u, material, quadrature_degree=2):
@@ -48,14 +49,14 @@ class MFrontNonlinearProblem(NonlinearProblem):
         self._init = True
 
         self.state_variables = []
-        self.gradients = {}
+        self.gradients = dict.fromkeys(self.material.get_gradient_names(), None)
         self.fluxes = {}
 
 
     def define_form(self):
         # residual form (internal forces)
-        self.L = sum([inner(g.variation(self.u_), f.function)
-                         for (f, g) in zip(self.fluxes.values(), self.gradients.values())])*self.axi_coeff*self.dx
+        self.L = sum([inner(g.variation(self.u_),f.function)
+                          for (f, g) in zip(self.fluxes.values(), self.gradients.values())])*self.axi_coeff*self.dx
         if self._Fext is not None:
             self.L -= self._Fext
         # tangent bilinear form
@@ -64,6 +65,8 @@ class MFrontNonlinearProblem(NonlinearProblem):
                          for (t, dg) in zip(f.tangent_blocks.values(), f.gradients)])*self.axi_coeff*self.dx
 
     def register_gradient(self, name, expression, symmetric=None):
+        pos = self.material.get_gradient_names().index(name)
+        size = self.material.get_gradient_sizes()[pos]
         self.gradients.update({name: Gradient(self.u, expression, name, symmetric=symmetric)})
 
     def set_loading(self, Fext):

mfront_wrapper/utils.py

@@ -86,7 +86,7 @@ def gradient(g):
     return nonsymmetric_tensor_to_vector(g)
 
 def get_quadrature_element(cell, degree, dim=0):
-    if dim in [0, 1, (), (0,), (1,), (0, 0), (0, 1), (1, 0)]:
+    if dim in [0, 1, (), (0,), (1,), (0, 0), (0, 1), (1, 0), (1, 1)]:
         return FiniteElement("Quadrature", cell, degree=degree,
                             quad_scheme='default')
     elif type(dim) == int: