Commit 14788799 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Added thermoelasticity example

parent 925dfdec
#!/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
length = 30e-3
width = 5.4e-3
mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)
V = FunctionSpace(mesh, "CG", 2)
T = Function(V, name="Temperature")
def left(x, on_boundary):
return near(x[1], 0) and on_boundary
def right(x, on_boundary):
return near(x[0], length) and on_boundary
Tl = 300
Tr = 800
bc = [DirichletBC(V, Constant(Tl), left),
DirichletBC(V, Constant(Tr), right)]
facets = MeshFunction("size_t", mesh, dim-1)
ds = Measure("ds", subdomain_data=facets)
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"StationaryHeatTransfer",
hypothesis="plane_strain")
problem = mf.MFrontNonlinearProblem(T, material)
problem.bc = bc
j = Function(Wjs, name="Current heat flux")
gT = Function(WgTs, name="Current temperature gradient")
dj_ddgT = Function(W_dj_ddgT_s, name="Derivative of the Heat Flux with respect to the temperature gradient")
dj_ddT = Function(W_dj_ddT_s, name="Derivative of the Heat Flux with respect to the temperature")
T0 = Function(V, name="Temperature at the beginning of the time step")
T1 = Function(V, name="Temperature estimate at the end of the end step")
dT = Function(V, name="Iteration correction")
v = TrialFunction(V)
T_ = TestFunction(V)
# definition of the stiffness matrix and residual using standard FEniCS operations
a_Newton = (inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))+
0*inner(grad(v), as_vector(dj_ddT))*T_)*dxm
# a_Newton = inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))*dxm
res = -inner(grad(T_), j)*dxm
problem.a = inner(self.strain_variation(self.du), dot(self.Ct, self.strain_variation(self.u_)))*measure*self.dx
problam.L = inner(self.strain_variation(self.u_), self.stress)*measure*self.dxl value of the deformation gradient
# using `MFront` conventions
T0.assign(Constant(Tl))
local_project(grad(T0), WgTs, gT)
# copy the strain values to `MGIS`
m.s0.gradients[:, :] = gT.vector().get_local().reshape((m.n, sj))
m.s1.gradients[:, :] = gT.vector().get_local().reshape((m.n, sj))
it = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
mgis_bv.integrate(m, it, 0, 0, m.n);
bs = sj*sj+sj*1
dj_ddgT.vector().set_local(np.lib.stride_tricks.as_strided(m.K,shape = (m.K.size // bs, sj, sj),
strides = (m.K.strides[1] * bs,
m.K.strides[1] * sj,
m.K.strides[1])).flatten())
dj_ddgT.vector().apply("insert")
dj_ddT.vector().set_local(np.lib.stride_tricks.as_strided(m.K[sj*sj:],
shape = (m.K.size // bs, sj),
strides = (m.K.strides[1] * bs,
m.K.strides[1])).flatten())
dj_ddT.vector().apply("insert")
print(dj_ddT.vector().norm("l2"))
T.assign(T0)
problem.solve(T.vector())
niter = 0
tol = 1.e-8
Nitermax = 200
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, dT.vector(), Res, "mumps")
# update the current estimate of the displacement at the end of the time step
T1.assign(T1+dT)
# compute the current estimate of the strain at the end of the
# time step using `MFront` conventions
local_project(grad(T1), WgTs, gT)
# copy the strain values to `MGIS`
m.s1.gradients[:, :] = gT.vector().get_local().reshape((m.n, sj))
# integrate the behaviour
it = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
# no time step here
mgis_bv.integrate(m, it, 0, 0, m.n);
# getting the stress and consistent tangent operator back to
# the FEniCS world.
j.vector().set_local(m.s1.thermodynamic_forces.flatten())
j.vector().apply("insert")
dj_ddgT.vector().set_local(np.lib.stride_tricks.as_strided(m.K,shape = (m.K.size // bs, sj, sj),
strides = (m.K.strides[1] * bs,
m.K.strides[1] * sj,
m.K.strides[1])).flatten())
dj_ddgT.vector().apply("insert")
dj_ddT.vector().set_local(np.lib.stride_tricks.as_strided(m.K[sj*sj:],
shape = (m.K.size // bs, sj),
strides = (m.K.strides[1] * bs,
m.K.strides[1])).flatten())
dj_ddT.vector().apply("insert")
# solve Newton system
A, Res = assemble_system(a_Newton, res, bc0)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
@DSL DefaultGenericBehaviour;
@Behaviour StationaryHeatTransfer;
@Author Thomas Helfer;
@Date 15/02/2019;
@Gradient TemperatureGradient gT;
gT.setGlossaryName("TemperatureGradient");
@Flux HeatFlux j;
j.setGlossaryName("HeatFlux");
@AdditionalTangentOperatorBlock dj_ddT;
@Parameter A = 0.0375;
@Parameter B = 2.165e-4;
@LocalVariable thermalconductivity k;
@Integrator{
const auto T_ = T + dT;
k = 1 / (A + B * T_);
j = k * (gT + dgT);
} // end of @Integrator
@TangentOperator {
dj_ddgT = k * tmatrix<N, N, real>::Id();
dj_ddT = -B * k * k * (gT + dgT);
} // end of @TangentOperator
/*!
* \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/StationaryHeatTransfer.hxx
* \brief this file implements the StationaryHeatTransfer Behaviour.
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 15 / 02 / 2019
*/
#ifndef LIB_TFELMATERIAL_STATIONARYHEATTRANSFER_HXX
#define LIB_TFELMATERIAL_STATIONARYHEATTRANSFER_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/Math/Vector/TVectorView.hxx"
#include"TFEL/Math/Matrix/TMatrixView.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/StationaryHeatTransferBehaviourData.hxx"
#include"TFEL/Material/StationaryHeatTransferIntegrationData.hxx"
#include "MFront/GenericBehaviour/State.hxx"
#include "MFront/GenericBehaviour/BehaviourData.hxx"
namespace tfel{
namespace material{
struct StationaryHeatTransferParametersInitializer
{
static StationaryHeatTransferParametersInitializer&
get();
double A;
double B;
double minimal_time_step_scaling_factor;
double maximal_time_step_scaling_factor;
void set(const char* const,const double);
/*!
* \brief convert a string to double
* \param[in] p : parameter
* \param[in] v : value
*/
static double getDouble(const std::string&,const std::string&);
private :
StationaryHeatTransferParametersInitializer();
StationaryHeatTransferParametersInitializer(const StationaryHeatTransferParametersInitializer&);
StationaryHeatTransferParametersInitializer&
operator=(const StationaryHeatTransferParametersInitializer&);
/*!
* \brief read the parameters from the given file
* \param[out] pi : parameters initializer
* \param[in] fn : file name
*/
static void readParameters(StationaryHeatTransferParametersInitializer&,const char* const);
};
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis,typename Type,bool use_qt>
class StationaryHeatTransfer;
//! \brief forward declaration
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream&,const StationaryHeatTransfer<hypothesis,Type,false>&);
/*!
* \class StationaryHeatTransfer
* \brief This class implements the StationaryHeatTransfer behaviour.
* \param hypothesis, modelling hypothesis.
* \param Type, numerical type.
* \author Thomas Helfer
* \date 15 / 02 / 2019
*/
template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
class StationaryHeatTransfer<hypothesis,Type,false> final
: public MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>,
public StationaryHeatTransferBehaviourData<hypothesis,Type,false>,
public StationaryHeatTransferIntegrationData<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 StationaryHeatTransfer&);
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 = tfel::math::tvector<(TVectorSize)*(TVectorSize)+(TVectorSize)*(1),real>;
using PhysicalConstants = tfel::PhysicalConstants<real>;
public :
typedef StationaryHeatTransferBehaviourData<hypothesis,Type,false> BehaviourData;
typedef StationaryHeatTransferIntegrationData<hypothesis,Type,false> IntegrationData;
typedef typename MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::SMFlag SMFlag;
typedef typename MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::SMType SMType;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::ELASTIC;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::SECANTOPERATOR;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::TANGENTOPERATOR;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::CONSISTENTTANGENTOPERATOR;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::NOSTIFFNESSREQUESTED;
using IntegrationResult = typename MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::IntegrationResult;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::SUCCESS;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::FAILURE;
using MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::UNRELIABLE_RESULTS;
private :
#line 17 "StationaryHeatTransfer.mfront"
thermalconductivity k;
#line 14 "StationaryHeatTransfer.mfront"
real A;
#line 15 "StationaryHeatTransfer.mfront"
real B;
real minimal_time_step_scaling_factor;
real maximal_time_step_scaling_factor;
//! Tangent operator;
TangentOperator Dt;
tfel::math::TMatrixView<N,N,real> dj_ddgT;
tfel::math::TVectorView<N,real> dj_ddT;
/*!
* \brief Update internal variables at end of integration
*/
void updateIntegrationVariables()
{}
/*!
* \brief Update internal variables at end of integration
*/
void updateStateVariables()
{}
/*!
* \brief Update auxiliary state variables at end of integration
*/
void updateAuxiliaryStateVariables()
{}
//! \brief Default constructor (disabled)
StationaryHeatTransfer() =delete ;
//! \brief Copy constructor (disabled)
StationaryHeatTransfer(const StationaryHeatTransfer&) = delete;
//! \brief Assignement operator (disabled)
StationaryHeatTransfer& operator = (const StationaryHeatTransfer&) = delete;
public:
/*!
* \brief Constructor
*/
StationaryHeatTransfer(const StationaryHeatTransferBehaviourData<hypothesis,Type,false>& src1,
const StationaryHeatTransferIntegrationData<hypothesis,Type,false>& src2)
: StationaryHeatTransferBehaviourData<hypothesis,Type,false>(src1),
StationaryHeatTransferIntegrationData<hypothesis,Type,false>(src2),
dj_ddgT(Dt.begin()),
dj_ddT(Dt.begin()+TVectorSize*TVectorSize)
{
using namespace std;
using namespace tfel::math;
using std::vector;
this->A = StationaryHeatTransferParametersInitializer::get().A;
this->B = StationaryHeatTransferParametersInitializer::get().B;
this->minimal_time_step_scaling_factor = StationaryHeatTransferParametersInitializer::get().minimal_time_step_scaling_factor;
this->maximal_time_step_scaling_factor = StationaryHeatTransferParametersInitializer::get().maximal_time_step_scaling_factor;
}
/*
* \brief constructor for the Generic interface
* \param[in] mgb_d: behaviour data
*/
StationaryHeatTransfer(const mfront::gb::BehaviourData& mgb_d)
: StationaryHeatTransferBehaviourData<hypothesis,Type,false>(mgb_d),
StationaryHeatTransferIntegrationData<hypothesis,Type,false>(mgb_d),
dj_ddgT(Dt.begin()),
dj_ddT(Dt.begin()+TVectorSize*TVectorSize)
{
using namespace std;
using namespace tfel::math;
using std::vector;
this->A = StationaryHeatTransferParametersInitializer::get().A;
this->B = StationaryHeatTransferParametersInitializer::get().B;
this->minimal_time_step_scaling_factor = StationaryHeatTransferParametersInitializer::get().minimal_time_step_scaling_factor;
this->maximal_time_step_scaling_factor = StationaryHeatTransferParametersInitializer::get().maximal_time_step_scaling_factor;
tfel::fsalgo::copy<TVectorSize >::exe(mgb_d.s0.gradients,this->gT.begin());
tfel::fsalgo::transform<TVectorSize>::exe(mgb_d.s1.gradients,mgb_d.s0.gradients,this->dgT.begin(),std::minus<real>());
tfel::fsalgo::copy<TVectorSize >::exe(mgb_d.s0.thermodynamic_forces,this->j.begin());
}
/*!
* \ brief initialize the behaviour with user code
*/
void initialize(){
using namespace std;
using namespace tfel::math;
using std::vector;
}
/*!
* \brief set the policy for "out of bounds" conditions
*/
void
setOutOfBoundsPolicy(const OutOfBoundsPolicy policy_value){
this->policy = policy_value;
} // end of setOutOfBoundsPolicy
/*!
* \return the modelling hypothesis
*/
constexpr ModellingHypothesis::Hypothesis
getModellingHypothesis() const{
return hypothesis;
} // end of getModellingHypothesis
/*!
* \brief check bounds
*/
void checkBounds() const{
} // end of checkBounds
IntegrationResult computePredictionOperator(const SMFlag,const SMType) override{
tfel::raise("StationaryHeatTransfer::computePredictionOperator: "
"unsupported prediction operator flag");
}
real getMinimalTimeStepScalingFactor() const override{
return this->minimal_time_step_scaling_factor;
}
std::pair<bool,real>
computeAPrioriTimeStepScalingFactor(const real current_time_step_scaling_factor) const override{
const auto time_scaling_factor = this->computeAPrioriTimeStepScalingFactorII();
return {time_scaling_factor.first,
std::min(std::min(std::max(time_scaling_factor.second,
this->minimal_time_step_scaling_factor),
this->maximal_time_step_scaling_factor),
current_time_step_scaling_factor)};
}
/*!
* \brief Integrate behaviour over the time step
*/
IntegrationResult
integrate(const SMFlag smflag, const SMType smt) override{
using namespace std;
using namespace tfel::math;
raise_if(smflag!=MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::STANDARDTANGENTOPERATOR,
"invalid tangent operator flag");
bool computeTangentOperator_ = smt!=NOSTIFFNESSREQUESTED;
#line 20 "StationaryHeatTransfer.mfront"
const auto T_ = this->T + this->dT;
#line 21 "StationaryHeatTransfer.mfront"
this->k = 1 / (this->A + this->B * T_);
#line 22 "StationaryHeatTransfer.mfront"
this->j = this->k * (this->gT + this->dgT);
this->updateIntegrationVariables();
this->updateStateVariables();
this->updateAuxiliaryStateVariables();
if(computeTangentOperator_){
if(!this->computeConsistentTangentOperator(smt)){
return MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::FAILURE;
}
}
return MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::SUCCESS;
}
std::pair<bool,real>
computeAPosterioriTimeStepScalingFactor(const real current_time_step_scaling_factor) const override{
const auto time_scaling_factor = this->computeAPosterioriTimeStepScalingFactorII();
return {time_scaling_factor.first,
std::min(std::min(std::max(time_scaling_factor.second,
this->minimal_time_step_scaling_factor),
this->maximal_time_step_scaling_factor),
current_time_step_scaling_factor)};
}
/*!
* \brief Update the internal energy at end of the time step
* \param[in] Psi_s: internal energy at end of the time step
*/
void computeInternalEnergy(real& Psi_s) const
{
Psi_s=0;
}
/*!
* \brief Update the dissipated energy at end of the time step
* \param[in] Psi_d: dissipated energy at end of the time step
*/
void computeDissipatedEnergy(real& Psi_d) const
{
Psi_d=0;
}
bool computeConsistentTangentOperator(const SMType smt){
using namespace std;
using namespace tfel::math;
using std::vector;
#line 26 "StationaryHeatTransfer.mfront"
this->dj_ddgT = this