Commit a6efe243 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Work on external state variables and non-linear heat transfer demo

parent ac4f20d8
#!/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", 1)
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, 1)
ds = Measure("ds", subdomain_data=facets)
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
"StationaryHeatTransfer",
hypothesis="plane_strain")
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
problem.bc = bc
# problem.register_gradient("TemperatureGradient", grad(T))
T.interpolate(Constant(Tl))
problem.solve(T.vector())
x = np.linspace(0, length)
import matplotlib.pyplot as plt
plt.plot(x, np.array([T(xi, width/2) for xi in x]))
#plt.figure()
#plot(project(flux.function[0], FunctionSpace(mesh, "DG",0)))
\ No newline at end of file
import matplotlib.pyplot as plt
from dolfin import *
import mfront_wrapper as mf
import numpy as np
......@@ -33,11 +34,10 @@ file_results.parameters["functions_share_mesh"] = True
selfweight = Expression(("0", "0", "-t*qmax"), t=0., qmax = 50e6, degree=0)
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"LogarithmicStrainPlasticity")
problem = mf.MFrontNonlinearProblem(u, material)
problem = mf.MFrontNonlinearProblem(u, material, bcs=bc)
problem.set_loading(dot(selfweight, u)*dx)
problem.bc = bc
p = problem.get_state_variable("EquivalentPlasticStrain")
assert (ufl.shape(p) == ())
......@@ -47,6 +47,7 @@ print(ufl.shape(epsel))
prm = problem.solver.parameters
prm["absolute_tolerance"] = 1e-6
prm["relative_tolerance"] = 1e-6
prm["linear_solver"] = "mumps"
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
......@@ -67,7 +68,6 @@ for (i, t) in enumerate(load_steps[1:]):
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")
......
......@@ -45,13 +45,12 @@ mat_prop = {"MatrixYoungModulus": 10.,
"FiberVolumeFraction": 0.01,
"Size": 1/16.}
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"MultiphaseModel",
hypothesis="plane_strain",
material_properties=mat_prop)
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
problem.bc = bc
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0, bcs=bc)
problem.register_gradient("MatrixStrain", sym(grad(u1)))
problem.register_gradient("FiberStrain", sym(grad(u2)))
problem.register_gradient("RelativeDisplacement", diag(u2-u1))
......
This diff is collapsed.
#!/usr/bin/env python
# coding: utf-8
# # Stationnary non-linear heat transfer
#
# ## Description of the non-linear constitutive heat transfer law
#
# ## FEniCS implementation
#
# We consider a rectanglar domain with imposed temperatures `Tl` (resp. `Tr`) on the left (resp. right boundaries). We want to solve for the temperature field `T` inside the domain using a $P^1$-interpolation. We initialize the temperature at value `Tl` throughout the domain. We finally load the material library with a `plane_strain` hypothesis.
# In[1]:
import matplotlib.pyplot as plt
from dolfin import *
import mfront_wrapper as mf
length = 30e-3
width = 5.4e-3
mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)
V = FunctionSpace(mesh, "CG", 1)
T = Function(V, name="Temperature")
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
Tl = 300
Tr = 800
T.interpolate(Constant(Tl))
bc = [DirichletBC(V, Constant(Tl), left),
DirichletBC(V, Constant(Tr), right)]
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"StationaryHeatTransfer",
hypothesis="plane_strain")
# The MFront behaviour implicitly declares the temperature as an external state variable called `"Temperature"`. We must therefore associate this external state variable to a known mechanical field. This can be achieved explicitly using the `register_external_state_variable` method. In the present case, this can be donc automatically since the name of the unknown temperature field matches the [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) name `"Temperature"`. In this case, the following message is printed:
# ```
# Automatic registration of 'Temperature' as an external state variable.
# ```
# For problems in which the temperature only acts as a parameter (no jacobian blocks with respect to the temperature), the temperature can be automatically registered as a constant value ($293.15 \text{ K}$ by default) or to any other (`dolfin.Constant` or `float`) value using the `register_external_state_variable` method.
#
# In the FEniCS interface, we instantiate the main mechanical unknown, here the temperature field `T` which has to be named `"Temperature"` in order to match MFront's predefined name. Using another name than this will later result in an error saying:
# ```
# ValueError: 'Temperature' could not be associated with a registered gradient or a known state variable.
# ```
#
# The MFront behaviour declares the field `"TemperatureGradient"` as a Gradient variable, with its associated Flux called `"HeatFlux"`. We can check that the `material` object retrieves MFront's gradient and flux names, as well as the different tangent operator blocks which have been defined, namely `dj_ddgT` and `dj_ddT` in the present case:
# In[2]:
print(material.get_gradient_names())
print(material.get_flux_names())
print(["d{}_d{}".format(*t) for t in material.get_tangent_block_names()])
# When defining the non-linear problem, we will specify the boundary conditions and the requested quadrature degree which will control the number of quadrature points used in each cell to compute the non-linear constitutive law. Here, we specify a quadrature of degree 2 (i.e. 3 Gauss points for a triangular element). Finally, we need to associate to MFront gradient object the corresponding UFL expression as a function of the unknown field `T`. To do so, we use the `register_gradient` method linking MFront `"TemperatureGradient"` object to the UFL expression `grad(T)`. Doing so, the corresponding non-linear variational problem will be automatically be built:
#
# \begin{equation}
# F(\widehat{T}) = \int_\Omega \boldsymbol{j}\cdot \nabla \widehat{T} \text{dx} = 0 \quad \forall \widehat{T}
# \end{equation}
# In[3]:
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=2, bcs=bc)
problem.register_gradient("TemperatureGradient", grad(T))
# From the two tangent operator blocks `dj_ddgT` and `dj_ddT`, it will automatically be deduced that the heat flux $\boldsymbol{j}$ is a function of both the temperature gradient $\boldsymbol{g}=\nabla T$ and the temperature itself i.e. $\boldsymbol{j}=\boldsymbol{j}(\boldsymbol{g}, T)$. The following tangent bilinear form will therefore be used when solving the above non-linear problem:
#
# \begin{equation}
# J(\widehat{T},T^*) = \int_{\Omega} \nabla \widehat{T}\cdot\left(\dfrac{\partial \boldsymbol{j}}{\partial \boldsymbol{g}}\cdot \nabla T^*+\dfrac{\partial \boldsymbol{j}}{\partial T}\cdot T^*\right) \text{dx}
# \end{equation}
#
# Similarly to the case of external state variables, common gradient expressions for some [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) names have been already predefined which avoid calling explicitly the `register_gradient` method. Predefined expressions can be obtained from:
# In[4]:
mf.list_predefined_gradients()
# We can see that the name `"Temperature Gradient"` is in fact a predefined gradient. Omitting calling the `register_gradient` method will in this case print the following message upon calling `solve`:
# ```
# Automatic registration of 'TemperatureGradient' as grad(Temperature).
# ```
# meaning that a predefined gradient name has been found and registered as the UFL expression $\nabla T$.
#
# We finally solve the non-linear problem using a default Newton non-linear solver. The `solve` method returns the number of Newton iterations (4 in the present case) and converged status .
# In[5]:
problem.solve(T.vector())
# We finally check that the thermal conductivity coefficient $k$, computed from the ratio between the horizontal heat flux and temperature gradient matches the temperature-dependent expressions implemented in the MFront behaviour.
# In[7]:
j = problem.fluxes["HeatFlux"].function
g = problem.gradients["TemperatureGradient"].function
k_gauss = j.vector().get_local()[::2]/g.vector().get_local()[::2]
T_gauss = problem.state_variables["external"]["Temperature"].function.vector().get_local()
A = 0.0375;
B = 2.165e-4;
k_ref = 1/(A + B*T_gauss)
plt.plot(T_gauss, k_gauss, 'o', label="FE")
plt.plot(T_gauss, k_ref, '.', label="ref")
plt.xlabel(r"Temperature $T\: (K)$")
plt.ylabel(r"Thermal conductivity $k\: (W.m^{-1}.K^{-1})$")
plt.legend()
plt.show()
# In[ ]:
......@@ -36,13 +36,12 @@ d = Function(Vd, name="Damage")
dold = Function(Vd, name="Previous damage")
material_u = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
material_u = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"PhaseFieldDisplacement",
hypothesis="plane_strain",
material_properties={"YoungModulus": 200.,
"PoissonRatio": 0.2})
problem_u = mf.MFrontNonlinearProblem(u, material_u, quadrature_degree=0)
problem_u.bc = bcu
problem_u = mf.MFrontNonlinearProblem(u, material_u, quadrature_degree=0, bcs=bcu)
problem_u.integration_type = mgis_bv.IntegrationType.IntegrationWithTangentOperator
d0 = problem_u.get_state_variable("Damage")
psi =problem_u.get_state_variable("DamageDrivingForce")
......@@ -52,18 +51,18 @@ prm_u["report"] = False
def solve_u(t):
Uimp.t = t
d0.interpolate(d)
problem_u.affect_state_variables()
problem_u.affect_internal_state_variables()
problem_u.solve(u.vector())
Gc = 1.
l0 = 0.02
material_d = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
material_d = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"PhaseFieldDamage",
hypothesis="plane_strain",
material_properties={"RegularizationLength": l0,
"FractureEnergy": Gc})
problem_d = mf.MFrontNonlinearProblem(d, material_d, quadrature_degree=0)
problem_d = mf.MFrontNonlinearProblem(d, material_d, quadrature_degree=0, bcs=bcd)
problem_d.register_gradient("Damage", d)
problem_d.register_gradient("DamageGradient", grad(d))
H = problem_d.get_state_variable("HistoryFunction")
......@@ -79,8 +78,7 @@ elif model == "with_threshold":
newH = ppos(psi-Gc/l0)
def solve_d(bc=[]):
H.assign(local_project(Max(H, newH), H.function_space(), problem_d.dx))
problem_d.affect_state_variables()
problem_d.bc = bcd
problem_d.affect_internal_state_variables()
dold.assign(d)
problem_d.solve(d.vector())
......
import matplotlib.pyplot as plt
from dolfin import *
import mfront_wrapper as mf
import numpy as np
......@@ -13,9 +14,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,12 +41,11 @@ 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 = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4, bcs=bc)
problem.set_loading(loading*dot(n, u)*measure*ds(4))
p = problem.get_state_variable("EquivalentPlasticStrain")
......@@ -59,7 +59,7 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic_strain")
Nincr = 40
Nincr = 10
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):
......@@ -72,8 +72,6 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
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}$")
......
@DSL Default;
@DSL DefaultDSL;
@Author Tran Thang DANG;
@Date 05/02/2016;
@Behaviour PhaseFieldDisplacement;
......
......@@ -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 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
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 PhaseFieldDamage-generic.cxx PhaseFieldDamage.cxx PhaseFieldDisplacement-generic.cxx PhaseFieldDisplacement.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 MultiphaseModel3D-generic.o MultiphaseModel3D.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 PhaseFieldDisplacement-generic.o PhaseFieldDisplacement.o PhaseFieldDamage-generic.o PhaseFieldDamage.o
@$(CXX) -shared $^ -o $@ -L"$(strip $(shell tfel-config --library-path))" $(patsubst %,-l%,$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants))
install :
......
......@@ -21,7 +21,11 @@ sources : {
"TensionCompressionSplit-generic.cxx",
"TensionCompressionSplit.cxx",
"MultiphaseModel3D-generic.cxx",
"MultiphaseModel3D.cxx"
"MultiphaseModel3D.cxx",
"PhaseFieldDisplacement-generic.cxx",
"PhaseFieldDisplacement.cxx",
"PhaseFieldDamage-generic.cxx",
"PhaseFieldDamage.cxx"
};
cppflags : {
"$(shell tfel-config --cppflags --compiler-flags)"
......@@ -68,7 +72,17 @@ epts : {
"TensionCompressionSplit_PlaneStrain",
"TensionCompressionSplit_GeneralisedPlaneStrain",
"TensionCompressionSplit_Tridimensional",
"MultiphaseModel3D_Tridimensional"
"MultiphaseModel3D_Tridimensional",
"PhaseFieldDisplacement_AxisymmetricalGeneralisedPlaneStrain",
"PhaseFieldDisplacement_Axisymmetrical",
"PhaseFieldDisplacement_PlaneStrain",
"PhaseFieldDisplacement_GeneralisedPlaneStrain",
"PhaseFieldDisplacement_Tridimensional",
"PhaseFieldDamage_AxisymmetricalGeneralisedPlaneStrain",
"PhaseFieldDamage_Axisymmetrical",
"PhaseFieldDamage_PlaneStrain",
"PhaseFieldDamage_GeneralisedPlaneStrain",
"PhaseFieldDamage_Tridimensional"
};
};
headers : {
......@@ -103,6 +117,14 @@ headers : {
"MFront/GenericBehaviour/MultiphaseModel3D-generic.hxx",
"TFEL/Material/MultiphaseModel3D.hxx",
"TFEL/Material/MultiphaseModel3DBehaviourData.hxx",
"TFEL/Material/MultiphaseModel3DIntegrationData.hxx"
"TFEL/Material/MultiphaseModel3DIntegrationData.hxx",
"MFront/GenericBehaviour/PhaseFieldDisplacement-generic.hxx",
"TFEL/Material/PhaseFieldDisplacement.hxx",
"TFEL/Material/PhaseFieldDisplacementBehaviourData.hxx",
"TFEL/Material/PhaseFieldDisplacementIntegrationData.hxx",
"MFront/GenericBehaviour/PhaseFieldDamage-generic.hxx",
"TFEL/Material/PhaseFieldDamage.hxx",
"TFEL/Material/PhaseFieldDamageBehaviourData.hxx",
"TFEL/Material/PhaseFieldDamageIntegrationData.hxx"
};
};
......@@ -9,8 +9,8 @@ from dolfin import *
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)
from .nonlinear_material import MFrontNonlinearMaterial
from .nonlinear_problem import MFrontNonlinearProblem
from .nonlinear_problem import MFrontNonlinearProblem, list_predefined_gradients
from .gradient_flux import *
\ No newline at end of file
......@@ -11,7 +11,7 @@ import ufl
from mfront_wrapper.utils import *
class Gradient:
def __init__(self, variable, expression, name=None, symmetric=None):
def __init__(self, variable, expression, name, symmetric=None):
self.variable = variable
if symmetric is None:
self.expression = expression
......@@ -33,10 +33,8 @@ class Gradient:
#self.expression = as_vector([self.expression])
else:
self.shape = shape
if name is None:
self.name = "Gradient-{}".format(np.random.randint(1e6))
else:
self.name = name
self.name = name
def __call__(self, v):
return ufl.replace(self.expression, {self.variable:v})
def variation(self, v):
......@@ -76,9 +74,26 @@ class Gradient:
# symmetric_tensor_to_vector(grad(variable) +
# Identity(variable.geometric_dimension()), T22=1), name)
class Var(Gradient):
""" The variable itself : u """
def __init__(self, variable, name=None):
Gradient.__init__(self, variable, variable, name)
def __init__(self, variable, name):
return Gradient.__init__(self, variable, variable, name)
# class Var:
# """ A variable """
# def __init__(self, expression, name=None):
# self.expression = expression
# shape = ufl.shape(self.expression)
# if len(shape)==1:
# self.shape = shape[0]
# elif shape==():
# self.shape = 0
# else:
# self.shape = shape
# if name is None:
# self.name = expression.name()
# else:
# self.name = name
# def variation(self, v):
# """ Directional derivative in direction v """
# return v
class Flux:
......
import mgis.behaviour as mgis_bv
from .gradient_flux import Var
import dolfin
mgis_hypothesis = {"plane_strain": mgis_bv.Hypothesis.PlaneStrain,
......@@ -9,15 +10,11 @@ mgis_hypothesis = {"plane_strain": mgis_bv.Hypothesis.PlaneStrain,
class MFrontNonlinearMaterial:
def __init__(self, path, name, hypothesis="3d",
material_properties={}, external_state_variables={"Temperature": 293.15}):
material_properties={}, external_state_variables={}):
self.path = path
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
self.material_properties = material_properties
self.external_state_variables = external_state_variables
# finite strain options
......@@ -25,16 +22,18 @@ class MFrontNonlinearMaterial:
bopts.stress_measure = mgis_bv.FiniteStrainBehaviourOptionsStressMeasure.PK1
bopts.tangent_operator = mgis_bv.FiniteStrainBehaviourOptionsTangentOperator.DPK1_DF
# Loading the behaviour
self.finite_strain = mgis_bv.isStandardFiniteStrainBehaviour(self.path, self.name)
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)
else:
self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis)
def set_data_manager(self, ngauss):
# Setting the material data manager
self.data_manager = mgis_bv.MaterialDataManager(self.behaviour, ngauss)
self.update_material_properties()
self.update_external_state_variable()
# self.update_external_state_variables()
def update_material_properties(self, material_properties=None):
if material_properties is not None:
......@@ -48,21 +47,33 @@ class MFrontNonlinearMaterial:
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:
self.external_state_variables = external_state_variables
def update_external_state_variables(self, external_state_variables):
for s in [self.data_manager.s0, self.data_manager.s1]:
for (key, value) in self.external_state_variables.items():
mgis_bv.setExternalStateVariable(s, key, value)
for (key, value) in external_state_variables.items():
if type(value) in [int, float]:
mgis_bv.setExternalStateVariable(s, key, value)
elif isinstance(value, dolfin.Constant):
mgis_bv.setExternalStateVariable(s, key, float(value))
else:
if isinstance(value, dolfin.Function):
value = value.vector().get_local()
elif isinstance(value, Var):
value.function.interpolate(value.expression)
value = value.function.vector().get_local()
mgis_bv.setExternalStateVariable(s, key, value, mgis_bv.MaterialStateManagerStorageMode.LocalStorage)
def get_state_variable_names(self):
def get_external_state_variable_names(self):
return [svar.name for svar in self.behaviour.external_state_variables]
def get_internal_state_variable_names(self):
return [svar.name for svar in self.behaviour.internal_state_variables]
def get_gradient_names(self):
return [svar.name for svar in self.behaviour.gradients]
def get_flux_names(self):
return [svar.name for svar in self.behaviour.thermodynamic_forces]
def get_state_variable_sizes(self):
def get_external_state_variable_sizes(self):
return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.external_state_variables]
def get_internal_state_variable_sizes(self):
return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.internal_state_variables]
def get_gradient_sizes(self):
return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.gradients]
......
......@@ -18,9 +18,17 @@ predefined_gradients = {
"TemperatureGradient": Grad,
"DeformationGradient": Id_Grad}
predefined_external_state_variables = {"Temperature": Constant(293.15)}
def list_predefined_gradients():
print("'TFEL gradient name' (Available hypotheses)")
print("---------------------------------------------")
for (g, value) in predefined_gradients.items():
print("{:22} {}".format("'{}'".format(g), tuple(str(v) for v in value.keys())))
print("")
class MFrontNonlinearProblem(NonlinearProblem):
def __init__(self, u, material, quadrature_degree=2):
def __init__(self, u, material, quadrature_degree=2, bcs=None):
NonlinearProblem.__init__(self)
self.u = u
self.V = self.u.function_space()
......@@ -32,7 +40,6 @@ class MFrontNonlinearProblem(NonlinearProblem):
self.integration_type = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
self.quadrature_degree = quadrature_degree
# self.set_quadrature_function_spaces()
cell = self.mesh.ufl_cell()
......@@ -40,7 +47,7 @@ class MFrontNonlinearProblem(NonlinearProblem):
# scalar quadrature space
self.W0 = FunctionSpace(self.mesh, W0e)
# compute Gauss points numbers
self.ngauss = self.W0.dim()
self.ngauss = Function(self.W0).vector().local_size()
# Set data manager
self.material.set_data_manager(self.ngauss)
......@@ -49,7 +56,10 @@ class MFrontNonlinearProblem(NonlinearProblem):
else:
assert self.u.geometric_dimension()==2, "Conflicting geometric dimension and material hypothesis"
self.bc = []
if bcs is None:
self.bcs = []
else:
self.bcs = bcs
self.dx = Measure("dx", metadata={"quadrature_degree": self.quadrature_degree,
"quadrature_scheme": "default"})
if self.axisymmetric:
......@@ -64,24 +74,34 @@ class MFrontNonlinearProblem(NonlinearProblem):
self._Fext = None
self._init = True
self.state_variables = dict.fromkeys(self.material.get_state_variable_names(), None)
self.state_variables = {"internal": dict.fromkeys(self.material.get_internal_state_variable_names(), None),
"external": dict.fromkeys(self.material.get_external_state_variable_names(), None)}
self.gradients = dict.fromkeys(self.material.get_gradient_names(), None)
self.fluxes = {}