Commit 88bca57f authored by Jeremy BLEYER's avatar Jeremy BLEYER

Working MFront wrapper for FEniCS

parents
from dolfin import *
import mfront_wrapper as mf
import numpy as np
length, width , height = 1., 40e-3, 100e-3
nb_elt_p, nb_elt_l = 10, 30
mesh = BoxMesh(Point(0, -width/2, -height/2.), Point(length, width/2, height/2.), nb_elt_l, 2, nb_elt_p)
V = VectorFunctionSpace(mesh, "CG", 2)
u = Function(V)
def left(x, on_boundary):
return near(x[0], 0) and on_boundary
def right(x, on_boundary):
return near(x[0], length) and on_boundary
bc = [DirichletBC(V, Constant((0.,)*3), left),
DirichletBC(V.sub(0), Constant(0.), right)]
facets = MeshFunction("size_t", mesh, 2)
AutoSubDomain(right).mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)
# Building function with vx=1 on left boundary to compute reaction
v = Function(V)
bcv = DirichletBC(V.sub(0), Constant(1.), left)
bcv.apply(v.vector())
Vpost = FunctionSpace(mesh, "CG", 1)
file_results = XDMFFile("beam_GL_plasticity.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
selfweight = Expression(("0", "0", "-t*qmax"), t=0., qmax = 50e6, degree=0)
material = mf.NonlinearMaterial('src/libBehaviour.so', "LogarithmicStrainPlasticity")
problem = mf.MFrontNonlinearProblem(u, material)
problem.set_loading(dot(selfweight, u)*dx)
problem.bc = bc
prm = problem.solver.parameters
prm["absolute_tolerance"] = 1e-6
prm["relative_tolerance"] = 1e-6
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
Nitermax, tol = 20, 1e-4 # parameters of the Newton-Raphson procedure
Nincr = 30
load_steps = np.linspace(0., 1., Nincr+1)
file_results.write(u, 0)
file_results.write(p_avg, 0)
results = np.zeros((Nincr+1, 3))
for (i, t) in enumerate(load_steps[1:]):
selfweight.t = t
problem.solve(u.vector())
results[i+1, 0] = assemble(-u[2]*ds(1))/(width*height)
results[i+1, 1] = t
results[i+1, 2] = assemble(action(-problem.L, v))
file_results.write(u, t)
# file_results.write(p_avg, t)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(results[:, 0], results[:, 1], "-o")
plt.xlabel(r"Displacement")
plt.ylabel("Load")
plt.figure()
plt.plot(results[:, 1], results[:, 2], "-o")
plt.xlabel(r"Load")
plt.ylabel("Horizontal Reaction")
plt.show()
from dolfin import *
import mfront_wrapper as mf
import numpy as np
hypothesis = "axisymmetric" # axisymmetric
Re, Ri = 1.3, 1. # external/internal radius
# elastic parameters
E = 70e3
nu = 0.3
# yield strength
sig0 = 250.
Et = E/100.
# hardening slope
H = E*Et/(E-Et)
mesh = Mesh("../meshes/thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "../meshes/thick_cylinder_facet_region.xml")
ds = Measure('ds', subdomain_data=facets)
if hypothesis == "plane_strain":
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)
measure = 1
elif hypothesis == "axisymmetric":
x = SpatialCoordinate(mesh)
q_lim = float(2*ln(Re/Ri)*sig0)
measure = 2*pi*abs(x[0])
V = VectorFunctionSpace(mesh, "CG", 2)
u = Function(V, name="Displacement")
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh)
loading = Expression("-q*t", q=q_lim, t=0, degree=2)
mat_prop = {"YoungModulus": E,
"PoissonRatio": nu,
"HardeningSlope": H,
"YieldStrength": sig0}
material = mf.NonlinearMaterial('materials/src/libBehaviour.so', 'IsotropicLinearHardeningPlasticity',
hypothesis=hypothesis,
material_properties=mat_prop)
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4)
problem.set_loading(loading*dot(n, u)*measure*ds(4))
problem.bc = bc
p = problem.register_state_variable(name="EquivalentPlasticStrain")
epsel = problem.register_state_variable(name="ElasticStrain", shape=4)
print(problem.state_variables)
file_results = XDMFFile("results/plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic_strain")
Nincr = 40
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
loading.t = t
problem.solve(u.vector())
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)
import matplotlib.pyplot as plt
plt.plot(results[:, 0], results[:, 1], "-o")
plt.xlabel("Displacement of inner boundary")
plt.ylabel(r"Applied pressure $q/q_{lim}$")
plt.show()
@DSL DefaultDSL;
@Behaviour IsotropicLinearHardeningPlasticity;
@Author Thomas Helfer;
@Date 14/10/2016;
@Description{
An implicit implementation of a simple
isotropic plasticity behaviour with
isotropic linear hardening.
The yield surface is defined by:
"\["
" f(\sigmaeq,p) = \sigmaeq-s_{0}-H\,p"
"\]"
}
@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");
@MaterialProperty stress H;
H.setEntryName("HardeningSlope");
@MaterialProperty stress s0;
s0.setGlossaryName("YieldStress");
@StateVariable StrainStensor eel;
eel.setGlossaryName("ElasticStrain");
@StateVariable strain p;
p.setGlossaryName("EquivalentPlasticStrain");
/*!
* computation of the prediction operator: we only provide the elastic
* operator.
*
* We could also provide a tangent operator, but this would mean
* saving an auxiliary state variable stating if a plastic loading
* occured at the previous time step.
*/
@PredictionOperator{
// silent "unused parameter" warning
static_cast<void>(smt);
const auto lambda = computeLambda(young,nu);
const auto mu = computeMu(young,nu);
Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
}
/*!
* behaviour integration using a fully implicit Euler-backwark scheme.
*/
@ProvidesSymmetricTangentOperator;
@Integrator{
const auto lambda = computeLambda(young,nu);
const auto mu = computeMu(young,nu);
eel += deto;
const auto se = 2*mu*deviator(eel);
const auto seq_e = sigmaeq(se);
const auto b = seq_e-s0-H*p>stress{0};
if(b){
const auto iseq_e = 1/seq_e;
const auto n = eval(3*se/(2*seq_e));
const auto cste = 1/(H+3*mu);
dp = (seq_e-s0-H*p)*cste;
eel -= dp*n;
if(computeTangentOperator_){
if(smt==CONSISTENTTANGENTOPERATOR){
Dt = (lambda*Stensor4::IxI()+2*mu*Stensor4::Id()
-4*mu*mu*(dp*iseq_e*(Stensor4::M()-(n^n))+cste*(n^n)));
} else {
Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
}
}
} else {
if(computeTangentOperator_){
Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
}
}
sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;
}
@DSL IsotropicPlasticMisesFlow;
@Behaviour LogarithmicStrainPlasticity;
@Author Helfer Thomas;
@Date 5/12/13;
@StrainMeasure Hencky;
@ElasticMaterialProperties {210e9,0.};
@Parameter s0 = 250e6;
@Parameter Et = 1e6;
@LocalVariable stress R;
@InitLocalVariables{
// hardening slope
R = young*Et/(young-Et);
}
@FlowRule{
f = seq-R*p-s0;
df_dseq = 1;
df_dp = -R;
}
/*!
* \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 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 */
/*!
* \file TFEL/Material/IsotropicLinearHardeningPlasticityBehaviourData.hxx
* \brief this file implements the IsotropicLinearHardeningPlasticityBehaviourData class.
* File generated by tfel version 3.3.0-dev
* \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/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 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;
}