diff --git a/demos/multiphase_model.py b/demos/multiphase_model.py index 75a8346951332be968ca2cff629ac7c25bd9d34b..61b2d4434e2d6f9c193c0b93147a386832e62225 100755 --- a/demos/multiphase_model.py +++ b/demos/multiphase_model.py @@ -10,9 +10,10 @@ import mfront_wrapper as mf import numpy as np from ufl import diag -width = 1. +width = 0.5 height = 1. -mesh = RectangleMesh(Point(0., 0.), Point(width, height), 10, 10) +mesh = RectangleMesh(Point(0., 0.), Point(width, height), 100, 100) +# mesh = BoxMesh(Point(0., 0., 0.), Point(width, height, height), 100, 2, 2) Ve = VectorElement("CG", mesh.ufl_cell(), 1) V = FunctionSpace(mesh, MixedElement([Ve, Ve])) @@ -37,13 +38,12 @@ bc = [DirichletBC(V.sub(0).sub(0), Constant(0), left), facets = MeshFunction("size_t", mesh, 1) ds = Measure("ds", subdomain_data=facets) - -mat_prop = {"MatrixYoungModulus": 70e3, - "MatrixPoissonRatio": 0.2, - "FiberYoungModulus": 200e3, +mat_prop = {"MatrixYoungModulus": 10., + "MatrixPoissonRatio": 0.45, + "FiberYoungModulus": 10000, "FiberPoissonRatio": 0.3, - "FiberVolumeFraction": 0.1, - "Size": 0.1} + "FiberVolumeFraction": 0.01, + "Size": 1/16.} material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", "MultiphaseModel", @@ -53,16 +53,40 @@ material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0) problem.bc = bc problem.register_gradient("MatrixStrain", sym(grad(u1))) -problem.register_gradient("InclusionStrain", sym(grad(u2))) +problem.register_gradient("FiberStrain", sym(grad(u2))) problem.register_gradient("RelativeDisplacement", diag(u2-u1)) +# E = Function(problem.W0) +# E.interpolate(Expression("20e3", degree=1)) +# mat_prop["MatrixYoungModulus"] = E +# problem.material.update_material_properties(mat_prop) problem.solve(u.vector()) -x = np.linspace(width/2, width, 21) +x = np.linspace(0, width, 21) import matplotlib.pyplot as plt 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"])) +plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "ob", label=r"$u_x^\textrm{matrix}$ EF") +plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "or", label=r"$u_x^\textrm{fiber}$ EF") +E1 = mat_prop["MatrixYoungModulus"] +nu1 = mat_prop["MatrixPoissonRatio"] +E2 = mat_prop["FiberYoungModulus"] +nu2 = mat_prop["FiberPoissonRatio"] +eta = mat_prop["FiberVolumeFraction"] +s = mat_prop["Size"] +mu1 = E1/2/(1+nu1) +mu2 = E2/2/(1+nu2) +lmbda1 = E1*nu1/(1+nu1)/(1-2*nu1) +lmbda2 = E2*nu2/(1+nu2)/(1-2*nu2) +kappa = 12*mu1*mu2/((1-eta)*mu2+eta*mu1)/s**2 +M1 = lmbda1+2*mu1 +M2 = lmbda2+2*mu2 +Er = eta*E2/(1-nu2**2) +ahom = M1+Er +l = sqrt(M1*Er/ahom/kappa) +um = lmbda1/ahom*(x+l*(Er/M1)*np.sinh(x/l)/np.cosh(width/l)) +ur = lmbda1/ahom*(x-l*np.sinh(x/l)/np.cosh(width/l)) +plt.plot(x, um, '--b', label=r"$u_x^\textrm{matrix}$ ref.") +plt.plot(x, ur, '--r', label=r"$u_x^\textrm{fiber}$ ref.") plt.legend() plt.show() \ No newline at end of file diff --git a/demos/phase_field_coupled.py b/demos/phase_field_coupled.py index 5b16660d8076463c1529813a43a92395bb747ced..aafc7d60e07874f50d1fa611ae6fbd6365f9ee5e 100755 --- a/demos/phase_field_coupled.py +++ b/demos/phase_field_coupled.py @@ -46,7 +46,6 @@ 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 diff --git a/materials/MultiphaseModel.mfront b/materials/MultiphaseModel.mfront deleted file mode 100644 index 4103f253071c07b502244d4091a25c03856b8250..0000000000000000000000000000000000000000 --- a/materials/MultiphaseModel.mfront +++ /dev/null @@ -1,73 +0,0 @@ -@DSL DefaultGenericBehaviour; -@Behaviour MultiphaseModel; -@ModellingHypotheses {PlaneStrain, Tridimensional}; - -@Gradient StrainStensor e₁; -e₁.setEntryName("MatrixStrain"); -@Flux StressStensor σ₁; -σ₁.setEntryName("MatrixStress"); - -@Gradient StrainStensor e₂; -e₂.setEntryName("FiberStrain"); -@Flux StressStensor σ₂; -σ₂.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₁; -@AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₂; -@AdditionalTangentOperatorBlock ∂I∕∂ΔV; - -@MaterialProperty stress Y1; -Y1.setEntryName("MatrixYoungModulus"); -@MaterialProperty real nu1; -nu1.setEntryName("MatrixPoissonRatio"); -@MaterialProperty stress Y2; -Y2.setEntryName("FiberYoungModulus"); -@MaterialProperty real nu2; -nu2.setEntryName("FiberPoissonRatio"); -@MaterialProperty real rho; -rho.setEntryName("FiberVolumeFraction"); -@MaterialProperty real s; -s.setEntryName("Size"); - -@ProvidesTangentOperator; -@Integrator { - // remove useless warnings, as we always compute the tangent operator - static_cast(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); - const auto μ₂ = rho*computeMu(Y2,nu2); - constexpr const auto λ₁₂ = 0.; - constexpr const auto μ₁₂ = 0.; - const auto kappa = μ₁/12/(1-rho)/pow(s,2); - 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₂ = C2; - σ₁ = ∂σ₁∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₁∕∂Δe₂ ⋅ (e₂ + Δe₂); - σ₂ = ∂σ₂∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₂∕∂Δe₂ ⋅ (e₂ + Δe₂); - I = kappa*(V + ΔV); - ∂I∕∂ΔV = kappa ⋅ I₄; -} diff --git a/materials/MultiphaseModel2D.mfront b/materials/MultiphaseModel2D.mfront new file mode 100644 index 0000000000000000000000000000000000000000..8ab2e67beeb3e5ef1d596d99507cefb2b40f081d --- /dev/null +++ b/materials/MultiphaseModel2D.mfront @@ -0,0 +1,83 @@ +@DSL DefaultGenericBehaviour; +@Behaviour MultiphaseModel; +@ModellingHypotheses {PlaneStrain}; + +@Gradient StrainStensor e₁; +e₁.setEntryName("MatrixStrain"); +@Flux StressStensor σ₁; +σ₁.setEntryName("MatrixStress"); + +@Gradient StrainStensor e₂; +e₂.setEntryName("FiberStrain"); +@Flux StressStensor σ₂; +σ₂.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₁; +@AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₂; +@AdditionalTangentOperatorBlock ∂I∕∂ΔV; + +@MaterialProperty stress Y1; +Y1.setEntryName("MatrixYoungModulus"); +@MaterialProperty real nu1; +nu1.setEntryName("MatrixPoissonRatio"); +@MaterialProperty stress Y2; +Y2.setEntryName("FiberYoungModulus"); +@MaterialProperty real nu2; +nu2.setEntryName("FiberPoissonRatio"); +@MaterialProperty real ρ; +ρ.setEntryName("FiberVolumeFraction"); +@MaterialProperty real s; +s.setEntryName("Size"); + +@ProvidesTangentOperator; +@Integrator { + // remove useless warnings, as we always compute the tangent operator + static_cast(computeTangentOperator_); + + const auto λ₁ = computeLambda(Y1,nu1); + const auto μ₁ = computeMu(Y1,nu1); + const auto λ₂ = computeLambda(Y2,nu2); + const auto μ₂ = computeMu(Y2,nu2); + const auto Eₒₑ¹ = λ₁+2*μ₁; + const auto Eₒₑ² = λ₂+2*μ₂; + const auto Eₒₑ = (1-ρ)*Eₒₑ¹ + ρ*Eₒₑ²; + const auto iEₒₑ = (1-ρ)/Eₒₑ¹ + ρ/Eₒₑ²; + const auto λ = (1-ρ)*λ₁ + ρ*λ₂; + const auto λiEₒₑ = (1-ρ)*λ₁/Eₒₑ¹ + ρ*λ₂/Eₒₑ²; + const auto λ2iEₒₑ = (1-ρ)*pow(λ₁,2)/Eₒₑ¹ + ρ*pow(λ₂,2)/Eₒₑ²; + const auto iμ = (1-ρ)/μ₁ + ρ/μ₂; + const auto κ = 12/iμ/pow(s,2); // interaction stiffness + const auto C₁₁₁₁ʰᵒᵐ = Eₒₑ - λ2iEₒₑ + pow(λiEₒₑ,2)/iEₒₑ; + const auto C₁₁₂₂ʰᵒᵐ = λiEₒₑ/iEₒₑ; + const Stensor4 Cʰᵒᵐ = { + C₁₁₁₁ʰᵒᵐ , C₁₁₂₂ʰᵒᵐ, C₁₁₂₂ʰᵒᵐ, 0., // + C₁₁₂₂ʰᵒᵐ, 1/iEₒₑ, λ, 0., // + C₁₁₂₂ʰᵒᵐ, λ, Eₒₑ, 0., // + 0., 0., 0., 2/iμ + }; + const auto C¹ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄; + const auto C² = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄; + const auto iΔC = invert(C²-C¹); + const auto Cᵃᵛᵍ = (1-ρ)⋅C¹ + ρ⋅C²; + const auto H = iΔC ⋅ (Cᵃᵛᵍ-Cʰᵒᵐ) ⋅ iΔC; + const auto D¹¹ = (1-ρ)⋅C¹-(C¹ ⋅ H ⋅ C¹); + const auto D²² = ρ⋅C²-(C² ⋅ H ⋅ C²); + const auto D¹² = C¹ ⋅ H ⋅ C²; + ∂σ₁∕∂Δe₁ = D¹¹; + ∂σ₁∕∂Δe₂ = D¹²; + ∂σ₂∕∂Δe₁ = transpose(D¹²); + ∂σ₂∕∂Δe₂ = D²²; + σ₁ = ∂σ₁∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₁∕∂Δe₂ ⋅ (e₂ + Δe₂); + σ₂ = ∂σ₂∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₂∕∂Δe₂ ⋅ (e₂ + Δe₂); + I = κ*(V + ΔV); + ∂I∕∂ΔV = κ ⋅ I₄; +} diff --git a/materials/MultiphaseModel3D.mfront b/materials/MultiphaseModel3D.mfront new file mode 100644 index 0000000000000000000000000000000000000000..e88c9be7343cea39b72b9a97d5b5e9213fe946f4 --- /dev/null +++ b/materials/MultiphaseModel3D.mfront @@ -0,0 +1,87 @@ +@DSL DefaultGenericBehaviour; +@Behaviour MultiphaseModel; +@ModellingHypotheses {Tridimensional}; + +@Gradient StrainStensor e₁; +e₁.setEntryName("MatrixStrain"); +@Flux StressStensor σ₁; +σ₁.setEntryName("MatrixStress"); + +@Gradient StrainStensor e₂; +e₂.setEntryName("FiberStrain"); +@Flux StressStensor σ₂; +σ₂.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₁; +@AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₂; +@AdditionalTangentOperatorBlock ∂I∕∂ΔV; + +@MaterialProperty stress Y1; +Y1.setEntryName("MatrixYoungModulus"); +@MaterialProperty real nu1; +nu1.setEntryName("MatrixPoissonRatio"); +@MaterialProperty stress Y2; +Y2.setEntryName("FiberYoungModulus"); +@MaterialProperty real nu2; +nu2.setEntryName("FiberPoissonRatio"); +@MaterialProperty real ρ; +ρ.setEntryName("FiberVolumeFraction"); +@MaterialProperty real s; +s.setEntryName("Size"); + +@ProvidesTangentOperator; +@Integrator { + // remove useless warnings, as we always compute the tangent operator + static_cast(computeTangentOperator_); + + const auto λ₁ = computeLambda(Y1,nu1); + const auto μ₁ = computeMu(Y1,nu1); + const auto λ₂ = computeLambda(Y2,nu2); + const auto μ₂ = computeMu(Y2,nu2); + const auto Eₒₑ¹ = λ₁+2*μ₁; + const auto Eₒₑ² = λ₂+2*μ₂; + const auto Eₒₑ = (1-ρ)*Eₒₑ¹ + ρ*Eₒₑ²; + const auto iEₒₑ = (1-ρ)/Eₒₑ¹ + ρ/Eₒₑ²; + const auto λ = (1-ρ)*λ₁ + ρ*λ₂; + const auto λiEₒₑ = (1-ρ)*λ₁/Eₒₑ¹ + ρ*λ₂/Eₒₑ²; + const auto λ2iEₒₑ = (1-ρ)*pow(λ₁,2)/Eₒₑ¹ + ρ*pow(λ₂,2)/Eₒₑ²; + const auto iμ = (1-ρ)/μ₁ + ρ/μ₂; + const auto μ = (1-ρ)*μ₁ + ρ*μ₂; + const auto κ = 12/iμ/pow(s,2); // interaction stiffness + const auto C₁₁₁₁ʰᵒᵐ = Eₒₑ - λ2iEₒₑ + pow(λiEₒₑ,2)/iEₒₑ; + const auto C₁₁₂₂ʰᵒᵐ = λiEₒₑ/iEₒₑ; + const Stensor4 Cʰᵒᵐ = { + C₁₁₁₁ʰᵒᵐ , C₁₁₂₂ʰᵒᵐ, C₁₁₂₂ʰᵒᵐ, 0., 0., 0., // + C₁₁₂₂ʰᵒᵐ, 1/iEₒₑ, λ, 0., 0., 0., // + C₁₁₂₂ʰᵒᵐ, λ, Eₒₑ, 0., 0., 0., // + 0., 0., 0., 2/iμ, 0., 0., // + 0., 0., 0., 0., 2/iμ, 0., // + 0., 0., 0., 0., 0., 2*μ // + }; + const auto C¹ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄; + const auto C² = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄; + const auto iΔC = invert(C²-C¹); + const auto Cᵃᵛᵍ = (1-ρ)*C¹ + ρ*C²; + const auto H = iΔC ⋅ (Cᵃᵛᵍ-Cʰᵒᵐ) ⋅ iΔC; + const auto D¹¹ = (1-ρ)*C¹-(C¹ ⋅ H ⋅ C¹); + const auto D²² = ρ*C²-(C² ⋅ H ⋅ C²); + const auto D¹² = C¹ ⋅ H ⋅ C²; + const auto D¹²ᵀ = C² ⋅ H ⋅ C¹; + ∂σ₁∕∂Δe₁ = D¹¹; + ∂σ₁∕∂Δe₂ = D¹²; + ∂σ₂∕∂Δe₁ = D¹²ᵀ; //transpose(D¹²); + ∂σ₂∕∂Δe₂ = D²²; + σ₁ = ∂σ₁∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₁∕∂Δe₂ ⋅ (e₂ + Δe₂); + σ₂ = ∂σ₂∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₂∕∂Δe₂ ⋅ (e₂ + Δe₂); + I = κ*(V + ΔV); + ∂I∕∂ΔV = κ ⋅ I₄; +} diff --git a/materials/src/Makefile.mfront b/materials/src/Makefile.mfront index a7825484c1b072e46a541c041025cf7041065933..87f2f5bd2d4c181b668037bae745f02d0c6bd3a9 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 TensionCompressionSplit-generic.cxx TensionCompressionSplit.cxx +SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiphaseModel-generic.cxx MultiphaseModel.cxx MultiphaseModel3D-generic.cxx MultiphaseModel3D.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 TensionCompressionSplit-generic.o TensionCompressionSplit.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 MultiphaseModel3D-generic.o MultiphaseModel3D.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 dbeb39ba6dd5346babd4dd365761e2cb79857056..da94076321df4672ed146548027d9da6bb3daa7a 100644 --- a/materials/src/targets.lst +++ b/materials/src/targets.lst @@ -19,7 +19,9 @@ sources : { "StationaryHeatTransfer-generic.cxx", "StationaryHeatTransfer.cxx", "TensionCompressionSplit-generic.cxx", -"TensionCompressionSplit.cxx" +"TensionCompressionSplit.cxx", +"MultiphaseModel3D-generic.cxx", +"MultiphaseModel3D.cxx" }; cppflags : { "$(shell tfel-config --cppflags --compiler-flags)" @@ -65,7 +67,8 @@ epts : { "TensionCompressionSplit_Axisymmetrical", "TensionCompressionSplit_PlaneStrain", "TensionCompressionSplit_GeneralisedPlaneStrain", -"TensionCompressionSplit_Tridimensional" +"TensionCompressionSplit_Tridimensional", +"MultiphaseModel3D_Tridimensional" }; }; headers : { @@ -96,6 +99,10 @@ headers : { "MFront/GenericBehaviour/TensionCompressionSplit-generic.hxx", "TFEL/Material/TensionCompressionSplit.hxx", "TFEL/Material/TensionCompressionSplitBehaviourData.hxx", -"TFEL/Material/TensionCompressionSplitIntegrationData.hxx" +"TFEL/Material/TensionCompressionSplitIntegrationData.hxx", +"MFront/GenericBehaviour/MultiphaseModel3D-generic.hxx", +"TFEL/Material/MultiphaseModel3D.hxx", +"TFEL/Material/MultiphaseModel3DBehaviourData.hxx", +"TFEL/Material/MultiphaseModel3DIntegrationData.hxx" }; }; diff --git a/mfront_wrapper/nonlinear_material.py b/mfront_wrapper/nonlinear_material.py index 4a55da32c2f125a38bf22480507007ed7a4bc964..f1804e17f9af01d7a9b8f92d0f5607220d6665f3 100644 --- a/mfront_wrapper/nonlinear_material.py +++ b/mfront_wrapper/nonlinear_material.py @@ -1,4 +1,5 @@ import mgis.behaviour as mgis_bv +import dolfin mgis_hypothesis = {"plane_strain": mgis_bv.Hypothesis.PlaneStrain, "plane_stress": mgis_bv.Hypothesis.PlaneStress, @@ -13,15 +14,21 @@ class MFrontNonlinearMaterial: self.name = name # Defining the modelling hypothesis self.hypothesis = mgis_hypothesis[hypothesis] - if self.hypothesis == mgis_bv.Hypothesis.Tridimensional: - self._tensorsize = 9 - else: - self._tensorsize = 5 + # if self.hypothesis == mgis_bv.Hypothesis.Tridimensional: + # self._tensorsize = 9 + # else: + # self._tensorsize = 5 self.material_properties = material_properties self.external_state_variables = external_state_variables + # finite strain options + bopts = mgis_bv.FiniteStrainBehaviourOptions() + bopts.stress_measure = mgis_bv.FiniteStrainBehaviourOptionsStressMeasure.PK1 + bopts.tangent_operator = mgis_bv.FiniteStrainBehaviourOptionsTangentOperator.DPK1_DF # Loading the behaviour self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis) self.finite_strain = self.behaviour.getBehaviourType()=="StandardFiniteStrainBehaviour" + if self.finite_strain: + self.behaviour = mgis_bv.load(bopts, self.path, self.name, self.hypothesis) def set_data_manager(self, ngauss): # Setting the material data manager @@ -34,7 +41,13 @@ class MFrontNonlinearMaterial: self.material_properties = material_properties for s in [self.data_manager.s0, self.data_manager.s1]: for (key, value) in self.material_properties.items(): - mgis_bv.setMaterialProperty(s, key, value) + if type(value) in [int, float]: + mgis_bv.setMaterialProperty(s, key, value) + else: + if isinstance(value, dolfin.Function): + value = value.vector().get_local() + mgis_bv.setMaterialProperty(s, key, value, mgis_bv.MaterialStateManagerStorageMode.LocalStorage) + def update_external_state_variable(self, external_state_variables=None): if external_state_variables is not None: @@ -54,12 +67,9 @@ class MFrontNonlinearMaterial: def get_gradient_sizes(self): return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.gradients] def get_flux_sizes(self): - return [self._tensorsize if (svar.name == "Stress" and self.finite_strain) else \ - mgis_bv.getVariableSize(svar, self.hypothesis) \ - for svar in self.behaviour.thermodynamic_forces] + return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.thermodynamic_forces] def get_tangent_block_names(self): return [(t[0].name, t[1].name) for t in self.behaviour.tangent_operator_blocks] def get_tangent_block_sizes(self): - return [tuple([self._tensorsize if (tt.name == "Stress" and self.finite_strain) \ - else mgis_bv.getVariableSize(tt, self.hypothesis) for tt in t]) \ + return [tuple([mgis_bv.getVariableSize(tt, self.hypothesis) for tt in t]) \ for t in self.behaviour.tangent_operator_blocks] \ No newline at end of file diff --git a/mfront_wrapper/nonlinear_problem.py b/mfront_wrapper/nonlinear_problem.py index eb6254b04428286dc1133c71ad5a26b9d5bb94f8..9a80df5e3dce71cd7976a24a4a4f0d122b356cb8 100644 --- a/mfront_wrapper/nonlinear_problem.py +++ b/mfront_wrapper/nonlinear_problem.py @@ -200,13 +200,13 @@ class MFrontNonlinearProblem(NonlinearProblem): f, g = block t = self.fluxes[f].tangent_blocks[g] block_shape = self.flattened_block_shapes[i] - if self.material.finite_strain and f == "Stress": - Ctv = t.vector().get_local() - mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager, - mgis_bv.FiniteStrainTangentOperator.DPK1_DF) - t.vector().set_local(Ctv) - else: - t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten()) + # if self.material.finite_strain and f == "Stress": + # Ctv = t.vector().get_local() + # # mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager, + # # mgis_bv.FiniteStrainTangentOperator.DPK1_DF) + # t.vector().set_local(Ctv) + # else: + t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten()) buff += block_shape def update_fluxes(self):