Commit ff873c92 authored by Jeremy BLEYER's avatar Jeremy BLEYER

More work on multiphase model

parent fb4e7fea
......@@ -12,7 +12,7 @@ from ufl import diag
width = 1.
height = 1.
mesh = RectangleMesh(Point(0., 0.), Point(width, height), 20, 10)
mesh = RectangleMesh(Point(0., 0.), Point(width, height), 10, 10)
Ve = VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([Ve, Ve]))
......@@ -38,9 +38,17 @@ facets = MeshFunction("size_t", mesh, 1)
ds = Measure("ds", subdomain_data=facets)
mat_prop = {"MatrixYoungModulus": 70e3,
"MatrixPoissonRatio": 0.2,
"FiberYoungModulus": 200e3,
"FiberPoissonRatio": 0.3,
"FiberVolumeFraction": 0.1,
"Size": 0.1}
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
"TwoPhasesGeneralizedElasticity",
hypothesis="plane_strain")
"MultiphaseModel",
hypothesis="plane_strain",
material_properties=mat_prop)
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
problem.bc = bc
......@@ -50,12 +58,24 @@ problem.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True)
problem.solve(u.vector())
u1 = u.sub(0, True)
u2 = u.sub(1, True)
x = np.linspace(0, width, 100)
x = np.linspace(width/2, width, 21)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x, np.array([u1(xi, height/2)[0] for xi in x]), label=r"$u_x$ (matrix)")
plt.plot(x, np.array([u2(xi, height/2)[0] for xi in x]), label=r"$u_x$ (inclusion)")
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
@DSL DefaultGenericBehaviour;
@Behaviour MultiPhase;
@Author Thomas Helfer;
@Date 15/02/2019;
@Gradient real eps1;
@Flux real sig1;
@Gradient real toto;
@Flux real sig2;
@Gradient real V;
@Flux real I;
// @AdditionalTangentOperatorBlock dsig1_deps1;
@Parameter C11 = 10e3;
@Parameter C12 = 5e-4;
@Parameter C22 = 2e3;
@Parameter kappa = 1e6;
@Integrator{
sig1 = C11*(eps1+deps1);
sig2 = C22*(toto+dtoto);
I = kappa*V
} // end of @Integrator
@TangentOperator {
//dj_ddgT = k * tmatrix<N, N, real>::Id();
dsig1_ddeps1 = C11* Stensor4::Id();
} // end of @TangentOperator
@DSL DefaultGenericBehaviour;
@Behaviour MultiphaseModel;
@ModellingHypotheses {PlaneStrain};
@Gradient StrainStensor e₁;
e₁.setEntryName("MatrixStrain");
@Flux StressStensor σ₁;
σ₁.setEntryName("MatrixStress");
@Gradient StrainStensor e₂;
e₂.setEntryName("InclusionStrain");
@Flux StressStensor σ₂;
σ₂.setEntryName("InclusionStress");
@Gradient StrainStensor V;
V.setEntryName("RelativeDisplacement");
@Flux StressStensor I;
I.setEntryName("InteractionForce");
@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<void>(computeTangentOperator_);
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);
∂σ₁∕∂Δe₁ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄;
∂σ₁∕∂Δe₂ = λ₁₂ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁₂ ⋅ I₄;
∂σ₂∕∂Δe₁ = ∂σ₁∕∂Δe₂;
∂σ₂∕∂Δe₂ = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄;
σ₁ = ∂σ₁∕∂Δe₁ ⋅ e₁ + ∂σ₁∕∂Δe₂ ⋅ e₂;
σ₂ = ∂σ₂∕∂Δe₁ ⋅ e₁ + ∂σ₂∕∂Δe₂ ⋅ e₂;
I = kappa*V;
∂I∕∂ΔV = kappa ⋅ I₄;
}
/*!
* \file IsotropicLinearHardeningPlasticity-generic.hxx
* \brief This file declares the umat interface for the IsotropicLinearHardeningPlasticity behaviour law
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
#ifndef LIB_GENERIC_ISOTROPICLINEARHARDENINGPLASTICITY_HXX
#define LIB_GENERIC_ISOTROPICLINEARHARDENINGPLASTICITY_HXX
#include"TFEL/Config/TFELConfig.hxx"
#include"MFront/GenericBehaviour/BehaviourData.h"
#ifdef _WIN32
#ifndef NOMINMAX
#define NOMINMAX
#endif /* NOMINMAX */
#include <windows.h>
#ifdef small
#undef small
#endif /* small */
#endif /* _WIN32 */
#ifndef MFRONT_SHAREDOBJ
#define MFRONT_SHAREDOBJ TFEL_VISIBILITY_EXPORT
#endif /* MFRONT_SHAREDOBJ */
#ifdef __cplusplus
extern "C"{
#endif /* __cplusplus */
MFRONT_SHAREDOBJ void
IsotropicLinearHardeningPlasticity_setOutOfBoundsPolicy(const int);
MFRONT_SHAREDOBJ int
IsotropicLinearHardeningPlasticity_setParameter(const char *const,const double);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int IsotropicLinearHardeningPlasticity_AxisymmetricalGeneralisedPlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int IsotropicLinearHardeningPlasticity_Axisymmetrical(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int IsotropicLinearHardeningPlasticity_PlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int IsotropicLinearHardeningPlasticity_GeneralisedPlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int IsotropicLinearHardeningPlasticity_Tridimensional(MFront_GB_BehaviourData* const);
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* LIB_GENERIC_ISOTROPICLINEARHARDENINGPLASTICITY_HXX */
/*!
* \file StationaryHeatTransfer-generic.hxx
* \brief This file declares the umat interface for the StationaryHeatTransfer behaviour law
* \author Thomas Helfer
* \date 15 / 02 / 2019
*/
#ifndef LIB_GENERIC_STATIONARYHEATTRANSFER_HXX
#define LIB_GENERIC_STATIONARYHEATTRANSFER_HXX
#include"TFEL/Config/TFELConfig.hxx"
#include"MFront/GenericBehaviour/BehaviourData.h"
#ifdef _WIN32
#ifndef NOMINMAX
#define NOMINMAX
#endif /* NOMINMAX */
#include <windows.h>
#ifdef small
#undef small
#endif /* small */
#endif /* _WIN32 */
#ifndef MFRONT_SHAREDOBJ
#define MFRONT_SHAREDOBJ TFEL_VISIBILITY_EXPORT
#endif /* MFRONT_SHAREDOBJ */
#ifdef __cplusplus
extern "C"{
#endif /* __cplusplus */
MFRONT_SHAREDOBJ void
StationaryHeatTransfer_setOutOfBoundsPolicy(const int);
MFRONT_SHAREDOBJ int
StationaryHeatTransfer_setParameter(const char *const,const double);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int StationaryHeatTransfer_AxisymmetricalGeneralisedPlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int StationaryHeatTransfer_Axisymmetrical(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int StationaryHeatTransfer_PlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int StationaryHeatTransfer_GeneralisedPlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int StationaryHeatTransfer_Tridimensional(MFront_GB_BehaviourData* const);
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* LIB_GENERIC_STATIONARYHEATTRANSFER_HXX */
/*!
* \file TFEL/Material/IsotropicLinearHardeningPlasticityBehaviourData.hxx
* \brief this file implements the IsotropicLinearHardeningPlasticityBehaviourData class.
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
#ifndef LIB_TFELMATERIAL_ISOTROPICLINEARHARDENINGPLASTICITY_BEHAVIOUR_DATA_HXX
#define LIB_TFELMATERIAL_ISOTROPICLINEARHARDENINGPLASTICITY_BEHAVIOUR_DATA_HXX
#include<limits>
#include<string>
#include<sstream>
#include<iostream>
#include<stdexcept>
#include<algorithm>
#include"TFEL/Raise.hxx"
#include"TFEL/PhysicalConstants.hxx"
#include"TFEL/Config/TFELConfig.hxx"
#include"TFEL/Config/TFELTypes.hxx"
#include"TFEL/Metaprogramming/StaticAssert.hxx"
#include"TFEL/TypeTraits/IsFundamentalNumericType.hxx"
#include"TFEL/TypeTraits/IsReal.hxx"
#include"TFEL/Math/General/IEEE754.hxx"
#include"TFEL/Math/stensor.hxx"
#include"TFEL/Math/Stensor/StensorConceptIO.hxx"
#include"TFEL/Math/tmatrix.hxx"
#include"TFEL/Math/Matrix/tmatrixIO.hxx"
#include"TFEL/Math/st2tost2.hxx"
#include"TFEL/Math/ST2toST2/ST2toST2ConceptIO.hxx"
#include"TFEL/Math/ST2toST2/ST2toST2View.hxx"
#include"TFEL/Material/ModellingHypothesis.hxx"
#include "MFront/GenericBehaviour/State.hxx"
#include "MFront/GenericBehaviour/BehaviourData.hxx"
namespace tfel{
namespace material{
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename,bool>
class IsotropicLinearHardeningPlasticityBehaviourData;
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename Type,bool use_qt>
class IsotropicLinearHardeningPlasticityIntegrationData;
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream&,const IsotropicLinearHardeningPlasticityBehaviourData<hypothesis,Type,false>&);
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
class IsotropicLinearHardeningPlasticityBehaviourData<hypothesis,Type,false>
{
static constexpr unsigned short N = ModellingHypothesisToSpaceDimension<hypothesis>::value;
TFEL_STATIC_ASSERT(N==1||N==2||N==3);
TFEL_STATIC_ASSERT(tfel::typetraits::IsFundamentalNumericType<Type>::cond);
TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<Type>::cond);
friend std::ostream& operator<< <>(std::ostream&,const IsotropicLinearHardeningPlasticityBehaviourData&);
/* integration data is declared friend to access driving variables at the beginning of the time step */
friend class IsotropicLinearHardeningPlasticityIntegrationData<hypothesis,Type,false>;
static constexpr unsigned short TVectorSize = N;
typedef tfel::math::StensorDimeToSize<N> StensorDimeToSize;
static constexpr unsigned short StensorSize = StensorDimeToSize::value;
typedef tfel::math::TensorDimeToSize<N> TensorDimeToSize;
static constexpr unsigned short TensorSize = TensorDimeToSize::value;
using ushort = unsigned short;
using Types = tfel::config::Types<N,Type,false>;
using real = typename Types::real;
using time = typename Types::time;
using length = typename Types::length;
using frequency = typename Types::frequency;
using stress = typename Types::stress;
using strain = typename Types::strain;
using strainrate = typename Types::strainrate;
using stressrate = typename Types::stressrate;
using temperature = typename Types::temperature;
using thermalexpansion = typename Types::thermalexpansion;
using thermalconductivity = typename Types::thermalconductivity;
using massdensity = typename Types::massdensity;
using energydensity = typename Types::energydensity;
using TVector = typename Types::TVector;
using Stensor = typename Types::Stensor;
using Stensor4 = typename Types::Stensor4;
using FrequencyStensor = typename Types::FrequencyStensor;
using ForceTVector = typename Types::ForceTVector;
using StressStensor = typename Types::StressStensor;
using StressRateStensor = typename Types::StressRateStensor;
using DisplacementTVector = typename Types::DisplacementTVector;
using StrainStensor = typename Types::StrainStensor;
using StrainRateStensor = typename Types::StrainRateStensor;
using StiffnessTensor = typename Types::StiffnessTensor;
using Tensor = typename Types::Tensor;
using FrequencyTensor = typename Types::FrequencyTensor;
using StressTensor = typename Types::StressTensor;
using ThermalExpansionCoefficientTensor = typename Types::ThermalExpansionCoefficientTensor;
using DeformationGradientTensor = typename Types::DeformationGradientTensor;
using DeformationGradientRateTensor = typename Types::DeformationGradientRateTensor;
using TemperatureGradient = typename Types::TemperatureGradient;
using HeatFlux = typename Types::HeatFlux;
using TangentOperator = StiffnessTensor;
using PhysicalConstants = tfel::PhysicalConstants<real>;
protected:
StrainStensor eto;
StressStensor sig;
#line 17 "IsotropicLinearHardeningPlasticity.mfront"
stress young;
#line 19 "IsotropicLinearHardeningPlasticity.mfront"
real nu;
#line 21 "IsotropicLinearHardeningPlasticity.mfront"
stress H;
#line 23 "IsotropicLinearHardeningPlasticity.mfront"
stress s0;
#line 26 "IsotropicLinearHardeningPlasticity.mfront"
StrainStensor eel;
#line 28 "IsotropicLinearHardeningPlasticity.mfront"
strain p;
temperature T;
public:
/*!
* \brief Default constructor
*/
IsotropicLinearHardeningPlasticityBehaviourData()
{}
/*!
* \brief copy constructor
*/
IsotropicLinearHardeningPlasticityBehaviourData(const IsotropicLinearHardeningPlasticityBehaviourData& src)
: eto(src.eto),
sig(src.sig),
young(src.young),
nu(src.nu),
H(src.H),
s0(src.s0),
eel(src.eel),
p(src.p),
T(src.T)
{}
/*
* \brief constructor for the Generic interface
* \param[in] mgb_d: behaviour data
*/
IsotropicLinearHardeningPlasticityBehaviourData(const mfront::gb::BehaviourData& mgb_d)
: young(mgb_d.s1.material_properties[0]),
nu(mgb_d.s1.material_properties[1]),
H(mgb_d.s1.material_properties[2]),
s0(mgb_d.s1.material_properties[3]),
eel(&mgb_d.s0.internal_state_variables[0]),
p(mgb_d.s0.internal_state_variables[StensorSize]),
T(mgb_d.s0.external_state_variables[0])
{
}
/*
* \brief Assignement operator
*/
IsotropicLinearHardeningPlasticityBehaviourData&
operator=(const IsotropicLinearHardeningPlasticityBehaviourData& src){
this->eto = src.eto;
this->sig = src.sig;
this->young = src.young;
this->nu = src.nu;
this->H = src.H;
this->s0 = src.s0;
this->eel = src.eel;
this->p = src.p;
this->T = src.T;
return *this;
}
void exportStateData(mfront::gb::State& mbg_s1) const
{
using namespace tfel::math;
tfel::fsalgo::copy<StensorSize>::exe(this->sig.begin(), mbg_s1.thermodynamic_forces);
tfel::fsalgo::copy<StensorSize>::exe(this->eel.begin(), mbg_s1.internal_state_variables);
mbg_s1.internal_state_variables[StensorSize] = this->p;
} // end of exportStateData
}; // end of IsotropicLinearHardeningPlasticityBehaviourDataclass
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const IsotropicLinearHardeningPlasticityBehaviourData<hypothesis,Type,false>& b)
{
os << "εᵗᵒ : " << b.eto << '\n';
os << "σ : " << b.sig << '\n';
os << "young : " << b.young << '\n';
os << "nu : " << b.nu << '\n';
os << "H : " << b.H << '\n';
os << "s0 : " << b.s0 << '\n';
os << "eel : " << b.eel << '\n';
os << "p : " << b.p << '\n';
os << "T : " << b.T << '\n';
return os;
}
} // end of namespace material
} // end of namespace tfel
#endif /* LIB_TFELMATERIAL_ISOTROPICLINEARHARDENINGPLASTICITY_BEHAVIOUR_DATA_HXX */
/*!
* \file TFEL/Material/IsotropicLinearHardeningPlasticityIntegrationData.hxx
* \brief this file implements the IsotropicLinearHardeningPlasticityIntegrationData class.
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
#ifndef LIB_TFELMATERIAL_ISOTROPICLINEARHARDENINGPLASTICITY_INTEGRATION_DATA_HXX
#define LIB_TFELMATERIAL_ISOTROPICLINEARHARDENINGPLASTICITY_INTEGRATION_DATA_HXX
#include<string>
#include<iostream>
#include<limits>
#include<stdexcept>
#include<algorithm>
#include"TFEL/Raise.hxx"
#include"TFEL/PhysicalConstants.hxx"
#include"TFEL/Config/TFELConfig.hxx"
#include"TFEL/Config/TFELTypes.hxx"
#include"TFEL/Metaprogramming/StaticAssert.hxx"
#include"TFEL/TypeTraits/IsFundamentalNumericType.hxx"
#include"TFEL/TypeTraits/IsScalar.hxx"
#include"TFEL/TypeTraits/IsReal.hxx"
#include"TFEL/TypeTraits/Promote.hxx"
#include"TFEL/Math/General/IEEE754.hxx"
#include"TFEL/Math/stensor.hxx"
#include"TFEL/Math/st2tost2.hxx"
#include "MFront/GenericBehaviour/State.hxx"
#include "MFront/GenericBehaviour/BehaviourData.hxx"
namespace tfel{
namespace material{
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename Type,bool use_qt>
class IsotropicLinearHardeningPlasticityIntegrationData;
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream&,const IsotropicLinearHardeningPlasticityIntegrationData<hypothesis,Type,false>&);
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
class IsotropicLinearHardeningPlasticityIntegrationData<hypothesis,Type,false>
{
static constexpr unsigned short N = ModellingHypothesisToSpaceDimension<hypothesis>::value;
TFEL_STATIC_ASSERT(N==1||N==2||N==3);
TFEL_STATIC_ASSERT(tfel::typetraits::IsFundamentalNumericType<Type>::cond);
TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<Type>::cond);
friend std::ostream& operator<< <>(std::ostream&,const IsotropicLinearHardeningPlasticityIntegrationData&);
static constexpr unsigned short TVectorSize = N;
typedef tfel::math::StensorDimeToSize<N> StensorDimeToSize;
static constexpr unsigned short StensorSize = StensorDimeToSize::value;
typedef tfel::math::TensorDimeToSize<N> TensorDimeToSize;
static constexpr unsigned short TensorSize = TensorDimeToSize::value;
using ushort = unsigned short;
using Types = tfel::config::Types<N,Type,false>;
using real = typename Types::real;
using time = typename Types::time;
using length = typename Types::length;
using frequency = typename Types::frequency;
using stress = typename Types::stress;
using strain = typename Types::strain;
using strainrate = typename Types::strainrate;
using stressrate = typename Types::stressrate;
using temperature = typename Types::temperature;
using thermalexpansion = typename Types::thermalexpansion;
using thermalconductivity = typename Types::thermalconductivity;