diff --git a/demos/phase_field_coupled.py b/demos/phase_field_coupled.py new file mode 100755 index 0000000000000000000000000000000000000000..5b16660d8076463c1529813a43a92395bb747ced --- /dev/null +++ b/demos/phase_field_coupled.py @@ -0,0 +1,98 @@ +#!/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 + + diff --git a/materials/MultiphaseModel.mfront b/materials/MultiphaseModel.mfront index 3d3931d5d2f763bc7f2a66bf710c305eeef01395..4103f253071c07b502244d4091a25c03856b8250 100644 --- a/materials/MultiphaseModel.mfront +++ b/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₄; } diff --git a/materials/TensionCompressionSplit.mfront b/materials/TensionCompressionSplit.mfront new file mode 100644 index 0000000000000000000000000000000000000000..503d1f11137e4009abbdc87731f0994103df404a --- /dev/null +++ b/materials/TensionCompressionSplit.mfront @@ -0,0 +1,71 @@ +@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; + } +} + diff --git a/materials/src/Makefile.mfront b/materials/src/Makefile.mfront index 5e46e5c67b65fffe58bc16e69ee19e2313307a54..a7825484c1b072e46a541c041025cf7041065933 100644 --- a/materials/src/Makefile.mfront +++ b/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 : diff --git a/materials/src/targets.lst b/materials/src/targets.lst index 93db986fa043811d7404f3dc9815a54b7685c600..dbeb39ba6dd5346babd4dd365761e2cb79857056 100644 --- a/materials/src/targets.lst +++ b/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" }; }; diff --git a/mfront_wrapper/gradient_flux.py b/mfront_wrapper/gradient_flux.py index 057c65c899c8f8b6f6371d9a945752aed68cfc71..68482f3294ce92c13af79ee1b215e313a3f46928 100644 --- a/mfront_wrapper/gradient_flux.py +++ b/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): diff --git a/mfront_wrapper/nonlinear_problem.py b/mfront_wrapper/nonlinear_problem.py index 9567d9d03bb759a265a2375cd661c4fbd2c7f204..eb6254b04428286dc1133c71ad5a26b9d5bb94f8 100644 --- a/mfront_wrapper/nonlinear_problem.py +++ b/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()):