Commit fb4e7fea authored by Jeremy BLEYER's avatar Jeremy BLEYER

Working tangent block formulation

parent 323faba3
......@@ -32,19 +32,13 @@ ds = Measure("ds", subdomain_data=facets)
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"HeatTransfer",
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
"StationaryHeatTransfer",
hypothesis="plane_strain")
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
problem.bc = bc
flux = mf.ThermalFlux(T)
problem.define_form(flux)
import mgis.behaviour as mgis_bv
for t in problem.material.behaviour.tangent_operator_blocks:
print('d{}_d{} size {}x{}'.format(t[0].name,t[1].name,
mgis_bv.getVariableSize(t[0], problem.material.hypothesis),
mgis_bv.getVariableSize(t[0], problem.material.hypothesis)))
problem.register_gradient("TemperatureGradient", grad(T))
T.interpolate(Constant(Tl))
......
#!/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
import numpy as np
from ufl import diag
width = 1.
height = 1.
mesh = RectangleMesh(Point(0., 0.), Point(width, height), 20, 10)
Ve = VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([Ve, Ve]))
u = Function(V, name="Displacements")
(u1, u2) = split(u)
def bottom(x, on_boundary):
return near(x[1], 0) and on_boundary
def top(x, on_boundary):
return near(x[1], height) and on_boundary
def left(x, on_boundary):
return near(x[0], 0) and on_boundary
bc = [DirichletBC(V.sub(0).sub(0), Constant(0), left),
DirichletBC(V.sub(1).sub(0), Constant(0), left),
DirichletBC(V.sub(0).sub(1), Constant(0), bottom),
DirichletBC(V.sub(1).sub(1), Constant(0), bottom),
DirichletBC(V.sub(0).sub(1), Constant(-1), top),
DirichletBC(V.sub(1).sub(1), Constant(-1), top)]
facets = MeshFunction("size_t", mesh, 1)
ds = Measure("ds", subdomain_data=facets)
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
"TwoPhasesGeneralizedElasticity",
hypothesis="plane_strain")
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
problem.bc = bc
problem.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True)
problem.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True)
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)
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.legend()
plt.show()
\ No newline at end of file
......@@ -13,9 +13,9 @@ Et = E/100.
# hardening slope
H = E*Et/(E-Et)
mesh = Mesh("../meshes/thick_cylinder.xml")
mesh = Mesh("meshes/thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "../meshes/thick_cylinder_facet_region.xml")
facets = MeshFunction("size_t", mesh, "meshes/thick_cylinder_facet_region.xml")
ds = Measure('ds', subdomain_data=facets)
......@@ -40,19 +40,20 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
"PoissonRatio": nu,
"HardeningSlope": H,
"YieldStrength": sig0}
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
"IsotropicLinearHardeningPlasticity",
hypothesis=hypothesis,
material_properties=mat_prop)
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4)
problem.bc = bc
problem.define_form(mf.Stress(u, name="sigma"))
# problem.define_form(mf.Stress(u, name="sigma"))
problem.register_gradient("Strain", sym(grad(u)), symmetric=True)
problem.set_loading(loading*dot(n, u)*measure*ds(4))
p = problem.get_state_variable(name="EquivalentPlasticStrain")
assert (ufl.shape(p) == ())
epsel = problem.get_state_variable(name="ElasticStrain")
assert (ufl.shape(epsel)==(4, ))
# p = problem.get_state_variable(name="EquivalentPlasticStrain")
# assert (ufl.shape(p) == ())
# epsel = problem.get_state_variable(name="ElasticStrain")
# assert (ufl.shape(epsel)==(4, ))
file_results = XDMFFile("results/plasticity_{}_results.xdmf".format(hypothesis))
file_results.parameters["flush_output"] = True
......@@ -69,9 +70,9 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
file_results.write(u, t)
p_avg.assign(project(epsel[0], P0))
file_results.write(p_avg, t)
results[i+1, :] = (u(Ri, 0)[0], t)
# p_avg.assign(project(epsel[0], P0))
# file_results.write(p_avg, t)
# results[i+1, :] = (u(Ri, 0)[0], t)
import matplotlib.pyplot as plt
......
@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
/*!
* \file LogarithmicStrainPlasticity-generic.hxx
* \brief This file declares the umat interface for the LogarithmicStrainPlasticity behaviour law
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
#ifndef LIB_GENERIC_LOGARITHMICSTRAINPLASTICITY_HXX
#define LIB_GENERIC_LOGARITHMICSTRAINPLASTICITY_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
LogarithmicStrainPlasticity_setOutOfBoundsPolicy(const int);
MFRONT_SHAREDOBJ int
LogarithmicStrainPlasticity_setParameter(const char *const,const double);
MFRONT_SHAREDOBJ int
LogarithmicStrainPlasticity_setUnsignedShortParameter(const char *const,const unsigned short);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int LogarithmicStrainPlasticity_AxisymmetricalGeneralisedPlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int LogarithmicStrainPlasticity_Axisymmetrical(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int LogarithmicStrainPlasticity_PlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int LogarithmicStrainPlasticity_GeneralisedPlaneStrain(MFront_GB_BehaviourData* const);
/*!
* \param[in,out] d: material data
*/
MFRONT_SHAREDOBJ int LogarithmicStrainPlasticity_Tridimensional(MFront_GB_BehaviourData* const);
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* LIB_GENERIC_LOGARITHMICSTRAINPLASTICITY_HXX */
......@@ -30,6 +30,7 @@
#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"
......@@ -138,12 +139,11 @@ IsotropicLinearHardeningPlasticityBehaviourData()
{}
/*!
* \brief Copy constructor
* \brief copy constructor
*/
IsotropicLinearHardeningPlasticityBehaviourData(const IsotropicLinearHardeningPlasticityBehaviourData& src)
: eto(src.eto),
sig(src.sig)
,
sig(src.sig),
young(src.young),
nu(src.nu),
H(src.H),
......
/*!
* \file TFEL/Material/LogarithmicStrainPlasticity.hxx
* \brief this file implements the LogarithmicStrainPlasticity Behaviour.
* File generated by tfel version 3.3.0
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
#ifndef LIB_TFELMATERIAL_LOGARITHMICSTRAINPLASTICITY_HXX
#define LIB_TFELMATERIAL_LOGARITHMICSTRAINPLASTICITY_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/IsReal.hxx"
#include"TFEL/Math/General/IEEE754.hxx"
#include"TFEL/Material/MaterialException.hxx"
#include"TFEL/Material/MechanicalBehaviour.hxx"
#include"TFEL/Material/MechanicalBehaviourTraits.hxx"
#include"TFEL/Material/OutOfBoundsPolicy.hxx"
#include"TFEL/Material/BoundsCheck.hxx"
#include"TFEL/Material/IsotropicPlasticity.hxx"
#include"TFEL/Material/Lame.hxx"
#include"TFEL/Material/Hosford1972YieldCriterion.hxx"
#include"TFEL/Material/LogarithmicStrainComputeAxialStrainIncrementElasticPrediction.hxx"
#include"TFEL/Material/LogarithmicStrainPlasticityBehaviourData.hxx"
#include"TFEL/Material/LogarithmicStrainPlasticityIntegrationData.hxx"
#include"TFEL/Math/General/BaseCast.hxx"
#include "MFront/GenericBehaviour/State.hxx"
#include "MFront/GenericBehaviour/BehaviourData.hxx"
namespace tfel{
namespace material{
struct LogarithmicStrainPlasticityParametersInitializer
{
static LogarithmicStrainPlasticityParametersInitializer&
get();
double young;
double nu;
double s0;
double Et;
double minimal_time_step_scaling_factor;
double maximal_time_step_scaling_factor;
double theta;
double epsilon;
unsigned short iterMax;
void set(const char* const,const double);
void set(const char* const,const unsigned short);
/*!
* \brief convert a string to double
* \param[in] p : parameter
* \param[in] v : value
*/
static double getDouble(const std::string&,const std::string&);
/*!
* \brief convert a string to unsigned short
* \param[in] p : parameter
* \param[in] v : value
*/
static unsigned short getUnsignedShort(const std::string&,const std::string&);
private :
LogarithmicStrainPlasticityParametersInitializer();
LogarithmicStrainPlasticityParametersInitializer(const LogarithmicStrainPlasticityParametersInitializer&);
LogarithmicStrainPlasticityParametersInitializer&
operator=(const LogarithmicStrainPlasticityParametersInitializer&);
/*!
* \brief read the parameters from the given file
* \param[out] pi : parameters initializer
* \param[in] fn : file name
*/
static void readParameters(LogarithmicStrainPlasticityParametersInitializer&,const char* const);
};
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis,typename Type,bool use_qt>
class LogarithmicStrainPlasticity;
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream&,const LogarithmicStrainPlasticity<hypothesis,Type,false>&);
/*!
* \class LogarithmicStrainPlasticity
* \brief This class implements the LogarithmicStrainPlasticity behaviour.
* \param hypothesis, modelling hypothesis.
* \param Type, numerical type.
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
class LogarithmicStrainPlasticity<hypothesis,Type,false> final
: public MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>,
public LogarithmicStrainPlasticityBehaviourData<hypothesis,Type,false>,
public LogarithmicStrainPlasticityIntegrationData<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 LogarithmicStrainPlasticity&);
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>;
public :
typedef LogarithmicStrainPlasticityBehaviourData<hypothesis,Type,false> BehaviourData;
typedef LogarithmicStrainPlasticityIntegrationData<hypothesis,Type,false> IntegrationData;
typedef typename MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::SMFlag SMFlag;
typedef typename MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::SMType SMType;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::ELASTIC;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::SECANTOPERATOR;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::TANGENTOPERATOR;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::CONSISTENTTANGENTOPERATOR;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::NOSTIFFNESSREQUESTED;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::STANDARDTANGENTOPERATOR;
using IntegrationResult = typename MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::IntegrationResult;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::SUCCESS;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::FAILURE;
using MechanicalBehaviour<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR,hypothesis,Type,false>::UNRELIABLE_RESULTS;
using StressFreeExpansionType = StrainStensor;
private :
typedef typename tfel::math::ComputeBinaryResult<strain,time,tfel::math::OpDiv>::Result DstrainDt;
typedef typename tfel::math::ComputeBinaryResult<DstrainDt,stress,tfel::math::OpDiv>::Result DF_DSEQ_TYPE;
StrainStensor deel;
strain dp;
temperature T_;
stress f;
real df_dseq;
stress df_dp;
StressStensor se;
stress seq;
stress seq_e;
StrainStensor n;
strain p_;
#line 13 "LogarithmicStrainPlasticity.mfront"
stress R;
stress lambda;
stress mu;
stress lambda_tdt;
stress mu_tdt;
stress young_tdt;
real nu_tdt;
real young;
real nu;
#line 10 "LogarithmicStrainPlasticity.mfront"
real s0;
#line 11 "LogarithmicStrainPlasticity.mfront"
real Et;
real minimal_time_step_scaling_factor;
real maximal_time_step_scaling_factor;
real theta;
real epsilon;
ushort iterMax;
//! Tangent operator;
TangentOperator Dt;
//! alias to the tangent operator;
TangentOperator& dsig_ddeto;
void computeFlow(){
using namespace std;
using namespace tfel::math;
using namespace tfel::material;
using std::vector;
#line 21 "LogarithmicStrainPlasticity.mfront"
this->f = this->seq-this->R*this->p_-this->s0;
#line 22 "LogarithmicStrainPlasticity.mfront"
this->df_dseq = 1;
#line 23 "LogarithmicStrainPlasticity.mfront"
this->df_dp = -this->R;
}
bool NewtonIntegration(){
using namespace std;
using namespace tfel::math;
bool converge=false;
bool inversible=true;
strain newton_f;
strain newton_df;
real newton_epsilon = 100*std::numeric_limits<real>::epsilon();
stress mu_3_theta = 3*(LogarithmicStrainPlasticity::theta)*(this->mu);
real surf;
unsigned int iter = 0u;
this->p_=this->p+this->dp;
while((converge==false)&&
(iter<this->iterMax)&&
(inversible==true)){
this->seq = std::max(this->seq_e-mu_3_theta*(this->dp),real(0.f));
this->computeFlow();
surf = (this->f)/(this->young);
if(((surf>newton_epsilon)&&((this->dp)>=0))||((this->dp)>newton_epsilon)){newton_f = surf;
newton_df = ((this->theta)*(this->df_dp)-mu_3_theta*(this->df_dseq))/(this->young);
} else {
newton_f =(this->dp);
newton_df = real(1.);
}
if(std::abs(base_cast(newton_df))>newton_epsilon){
this->dp -= newton_f/newton_df;
this->p_ = this->p + (this->theta)*(this->dp);
iter+=1;
converge = (std::abs(tfel::math::base_cast(newton_f))<this->epsilon);
} else {
inversible=false;
}
}
if(inversible==false){
return false;
}
if(iter==this->iterMax){
return false;
}
return true;
}
/*!
* \brief Update internal variables at end of integration
*/
void updateIntegrationVariables(){
}
/*!
* \brief Update internal variables at end of integration
*/
void updateStateVariables(){
this->eel += this->deel;
this->p += this->dp;
}
/*!
* \brief Update auxiliary state variables at end of integration
*/
void updateAuxiliaryStateVariables()
{}
//! \brief Default constructor (disabled)
LogarithmicStrainPlasticity() =delete ;
//! \brief Copy constructor (disabled)
LogarithmicStrainPlasticity(const LogarithmicStrainPlasticity&) = delete;
//! \brief Assignement operator (disabled)
LogarithmicStrainPlasticity& operator = (const LogarithmicStrainPlasticity&) = delete;
public: