Commit 12c69475 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Automatic registration and modifs on size

parent 83a2967f
...@@ -38,7 +38,7 @@ material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", ...@@ -38,7 +38,7 @@ material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0) problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
problem.bc = bc problem.bc = bc
problem.register_gradient("TemperatureGradient", grad(T)) # problem.register_gradient("TemperatureGradient", grad(T))
T.interpolate(Constant(Tl)) T.interpolate(Constant(Tl))
......
...@@ -33,15 +33,15 @@ file_results.parameters["functions_share_mesh"] = True ...@@ -33,15 +33,15 @@ file_results.parameters["functions_share_mesh"] = True
selfweight = Expression(("0", "0", "-t*qmax"), t=0., qmax = 50e6, degree=0) 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") "LogarithmicStrainPlasticity")
problem = mf.MFrontNonlinearProblem(u, material) problem = mf.MFrontNonlinearProblem(u, material)
problem.set_loading(dot(selfweight, u)*dx) problem.set_loading(dot(selfweight, u)*dx)
problem.bc = bc problem.bc = bc
p = problem.get_state_variable(name="EquivalentPlasticStrain") p = problem.get_state_variable("EquivalentPlasticStrain")
assert (ufl.shape(p) == ()) assert (ufl.shape(p) == ())
epsel = problem.get_state_variable(name="ElasticStrain") epsel = problem.get_state_variable("ElasticStrain")
print(ufl.shape(epsel)) print(ufl.shape(epsel))
prm = problem.solver.parameters prm = problem.solver.parameters
......
...@@ -52,9 +52,9 @@ material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", ...@@ -52,9 +52,9 @@ material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0) problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
problem.bc = bc problem.bc = bc
problem.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True) problem.register_gradient("MatrixStrain", sym(grad(u1)))
problem.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True) problem.register_gradient("InclusionStrain", sym(grad(u2)))
problem.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True) problem.register_gradient("RelativeDisplacement", diag(u2-u1))
problem.solve(u.vector()) problem.solve(u.vector())
......
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
...@@ -32,6 +33,7 @@ problem.bc = bc ...@@ -32,6 +33,7 @@ problem.bc = bc
problem.register_gradient("Damage", d) problem.register_gradient("Damage", d)
problem.register_gradient("DamageGradient", grad(d)) problem.register_gradient("DamageGradient", grad(d))
problem.solve(d.vector()) problem.solve(d.vector())
#TODO: manage state variables communication with external computation
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
p=plot(d) p=plot(d)
......
...@@ -46,14 +46,12 @@ for hypothesis in ["plane_strain", "axisymmetric"]: ...@@ -46,14 +46,12 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
material_properties=mat_prop) material_properties=mat_prop)
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4) problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4)
problem.bc = bc problem.bc = bc
# 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)) problem.set_loading(loading*dot(n, u)*measure*ds(4))
# p = problem.get_state_variable(name="EquivalentPlasticStrain") p = problem.get_state_variable("EquivalentPlasticStrain")
# assert (ufl.shape(p) == ()) assert (ufl.shape(p) == ())
# epsel = problem.get_state_variable(name="ElasticStrain") epsel = problem.get_state_variable("ElasticStrain")
# assert (ufl.shape(epsel)==(4, )) assert (ufl.shape(epsel)==(4, ))
file_results = XDMFFile("results/plasticity_{}_results.xdmf".format(hypothesis)) file_results = XDMFFile("results/plasticity_{}_results.xdmf".format(hypothesis))
file_results.parameters["flush_output"] = True file_results.parameters["flush_output"] = True
...@@ -70,9 +68,9 @@ for hypothesis in ["plane_strain", "axisymmetric"]: ...@@ -70,9 +68,9 @@ for hypothesis in ["plane_strain", "axisymmetric"]:
file_results.write(u, t) file_results.write(u, t)
# p_avg.assign(project(epsel[0], P0)) p_avg.assign(project(p, P0))
# file_results.write(p_avg, t) file_results.write(p_avg, t)
# results[i+1, :] = (u(Ri, 0)[0], t) results[i+1, :] = (u(Ri, 0)[0], t)
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
......
...@@ -12,7 +12,7 @@ INCLUDES := -I../include \ ...@@ -12,7 +12,7 @@ INCLUDES := -I../include \
CXXFLAGS := -Wall -Wfatal-errors -ansi $(shell tfel-config --oflags) -fPIC $(INCLUDES) CXXFLAGS := -Wall -Wfatal-errors -ansi $(shell tfel-config --oflags) -fPIC $(INCLUDES)
SRCCXX = MultiphaseModel-generic.cxx MultiphaseModel.cxx PhaseField-generic.cxx PhaseField.cxx SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiphaseModel-generic.cxx MultiphaseModel.cxx PhaseField-generic.cxx PhaseField.cxx PhaseFieldCoupled-generic.cxx PhaseFieldCoupled.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx
makefiles1 = $(SRCCXX:.cxx=.d) makefiles1 = $(SRCCXX:.cxx=.d)
makefiles2 = $(makefiles1:.cpp=.d) makefiles2 = $(makefiles1:.cpp=.d)
...@@ -22,7 +22,7 @@ makefiles = $(makefiles2) ...@@ -22,7 +22,7 @@ makefiles = $(makefiles2)
all : libBehaviour.so all : libBehaviour.so
libBehaviour.so : MultiphaseModel-generic.o MultiphaseModel.o PhaseField-generic.o PhaseField.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
@$(CXX) -shared $^ -o $@ -L"$(strip $(shell tfel-config --library-path))" $(patsubst %,-l%,$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants)) @$(CXX) -shared $^ -o $@ -L"$(strip $(shell tfel-config --library-path))" $(patsubst %,-l%,$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants))
install : install :
......
...@@ -9,7 +9,15 @@ sources : { ...@@ -9,7 +9,15 @@ sources : {
"MultiphaseModel-generic.cxx", "MultiphaseModel-generic.cxx",
"MultiphaseModel.cxx", "MultiphaseModel.cxx",
"PhaseField-generic.cxx", "PhaseField-generic.cxx",
"PhaseField.cxx" "PhaseField.cxx",
"PhaseFieldCoupled-generic.cxx",
"PhaseFieldCoupled.cxx",
"IsotropicLinearHardeningPlasticity-generic.cxx",
"IsotropicLinearHardeningPlasticity.cxx",
"LogarithmicStrainPlasticity-generic.cxx",
"LogarithmicStrainPlasticity.cxx",
"StationaryHeatTransfer-generic.cxx",
"StationaryHeatTransfer.cxx"
}; };
cppflags : { cppflags : {
"$(shell tfel-config --cppflags --compiler-flags)" "$(shell tfel-config --cppflags --compiler-flags)"
...@@ -30,7 +38,27 @@ epts : { ...@@ -30,7 +38,27 @@ epts : {
"PhaseField_Axisymmetrical", "PhaseField_Axisymmetrical",
"PhaseField_PlaneStrain", "PhaseField_PlaneStrain",
"PhaseField_GeneralisedPlaneStrain", "PhaseField_GeneralisedPlaneStrain",
"PhaseField_Tridimensional" "PhaseField_Tridimensional",
"PhaseFieldCoupled_AxisymmetricalGeneralisedPlaneStrain",
"PhaseFieldCoupled_Axisymmetrical",
"PhaseFieldCoupled_PlaneStrain",
"PhaseFieldCoupled_GeneralisedPlaneStrain",
"PhaseFieldCoupled_Tridimensional",
"IsotropicLinearHardeningPlasticity_AxisymmetricalGeneralisedPlaneStrain",
"IsotropicLinearHardeningPlasticity_Axisymmetrical",
"IsotropicLinearHardeningPlasticity_PlaneStrain",
"IsotropicLinearHardeningPlasticity_GeneralisedPlaneStrain",
"IsotropicLinearHardeningPlasticity_Tridimensional",
"LogarithmicStrainPlasticity_AxisymmetricalGeneralisedPlaneStrain",
"LogarithmicStrainPlasticity_Axisymmetrical",
"LogarithmicStrainPlasticity_PlaneStrain",
"LogarithmicStrainPlasticity_GeneralisedPlaneStrain",
"LogarithmicStrainPlasticity_Tridimensional",
"StationaryHeatTransfer_AxisymmetricalGeneralisedPlaneStrain",
"StationaryHeatTransfer_Axisymmetrical",
"StationaryHeatTransfer_PlaneStrain",
"StationaryHeatTransfer_GeneralisedPlaneStrain",
"StationaryHeatTransfer_Tridimensional"
}; };
}; };
headers : { headers : {
...@@ -41,6 +69,22 @@ headers : { ...@@ -41,6 +69,22 @@ headers : {
"MFront/GenericBehaviour/PhaseField-generic.hxx", "MFront/GenericBehaviour/PhaseField-generic.hxx",
"TFEL/Material/PhaseField.hxx", "TFEL/Material/PhaseField.hxx",
"TFEL/Material/PhaseFieldBehaviourData.hxx", "TFEL/Material/PhaseFieldBehaviourData.hxx",
"TFEL/Material/PhaseFieldIntegrationData.hxx" "TFEL/Material/PhaseFieldIntegrationData.hxx",
"MFront/GenericBehaviour/PhaseFieldCoupled-generic.hxx",
"TFEL/Material/PhaseFieldCoupled.hxx",
"TFEL/Material/PhaseFieldCoupledBehaviourData.hxx",
"TFEL/Material/PhaseFieldCoupledIntegrationData.hxx",
"MFront/GenericBehaviour/IsotropicLinearHardeningPlasticity-generic.hxx",
"TFEL/Material/IsotropicLinearHardeningPlasticity.hxx",
"TFEL/Material/IsotropicLinearHardeningPlasticityBehaviourData.hxx",
"TFEL/Material/IsotropicLinearHardeningPlasticityIntegrationData.hxx",
"MFront/GenericBehaviour/LogarithmicStrainPlasticity-generic.hxx",
"TFEL/Material/LogarithmicStrainPlasticity.hxx",
"TFEL/Material/LogarithmicStrainPlasticityBehaviourData.hxx",
"TFEL/Material/LogarithmicStrainPlasticityIntegrationData.hxx",
"MFront/GenericBehaviour/StationaryHeatTransfer-generic.hxx",
"TFEL/Material/StationaryHeatTransfer.hxx",
"TFEL/Material/StationaryHeatTransferBehaviourData.hxx",
"TFEL/Material/StationaryHeatTransferIntegrationData.hxx"
}; };
}; };
...@@ -16,14 +16,20 @@ class Gradient: ...@@ -16,14 +16,20 @@ class Gradient:
if symmetric is None: if symmetric is None:
self.expression = expression self.expression = expression
elif symmetric: elif symmetric:
if self.variable.geometric_dimension()==2:
self.expression = as_vector([symmetric_tensor_to_vector(expression)[i] for i in range(4)])
else:
self.expression = symmetric_tensor_to_vector(expression) self.expression = symmetric_tensor_to_vector(expression)
else:
if self.variable.geometric_dimension()==2:
self.expression = as_vector([nonsymmetric_tensor_to_vector(expression)[i] for i in range(5)])
else: else:
self.expression = nonsymmetric_tensor_to_vector(expression) self.expression = nonsymmetric_tensor_to_vector(expression)
shape = ufl.shape(self.expression) shape = ufl.shape(self.expression)
if len(shape)==1: if len(shape)==1:
self.shape = shape[0] self.shape = shape[0]
elif shape==(): elif shape==():
self.shape = 1 self.shape = 0
#self.expression = as_vector([self.expression]) #self.expression = as_vector([self.expression])
else: else:
self.shape = shape self.shape = shape
......
...@@ -13,10 +13,15 @@ class MFrontNonlinearMaterial: ...@@ -13,10 +13,15 @@ class MFrontNonlinearMaterial:
self.name = name self.name = name
# Defining the modelling hypothesis # Defining the modelling hypothesis
self.hypothesis = mgis_hypothesis[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.material_properties = material_properties
self.external_state_variables = external_state_variables self.external_state_variables = external_state_variables
# Loading the behaviour # Loading the behaviour
self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis) self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis)
self.finite_strain = self.behaviour.getBehaviourType()=="StandardFiniteStrainBehaviour"
def set_data_manager(self, ngauss): def set_data_manager(self, ngauss):
# Setting the material data manager # Setting the material data manager
...@@ -49,8 +54,12 @@ class MFrontNonlinearMaterial: ...@@ -49,8 +54,12 @@ class MFrontNonlinearMaterial:
def get_gradient_sizes(self): def get_gradient_sizes(self):
return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.gradients] return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.gradients]
def get_flux_sizes(self): def get_flux_sizes(self):
return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.thermodynamic_forces] return [self._tensorsize if (svar.name == "Stress" and self.finite_strain) else \
mgis_bv.getVariableSize(svar, self.hypothesis) \
for svar in self.behaviour.thermodynamic_forces]
def get_tangent_block_names(self): def get_tangent_block_names(self):
return [(t[0].name, t[1].name) for t in self.behaviour.tangent_operator_blocks] return [(t[0].name, t[1].name) for t in self.behaviour.tangent_operator_blocks]
def get_tangent_block_sizes(self): def get_tangent_block_sizes(self):
return [tuple([mgis_bv.getVariableSize(tt, self.hypothesis) for tt in t]) for t in self.behaviour.tangent_operator_blocks] return [tuple([self._tensorsize if (tt.name == "Stress" and self.finite_strain) \
\ No newline at end of file else mgis_bv.getVariableSize(tt, self.hypothesis) for tt in t]) \
for t in self.behaviour.tangent_operator_blocks]
\ No newline at end of file
...@@ -4,6 +4,21 @@ from mfront_wrapper.gradient_flux import * ...@@ -4,6 +4,21 @@ from mfront_wrapper.gradient_flux import *
import mgis.behaviour as mgis_bv import mgis.behaviour as mgis_bv
from ufl import shape from ufl import shape
SymGrad = {mgis_bv.Hypothesis.Tridimensional: lambda u: sym(grad(u)),
mgis_bv.Hypothesis.PlaneStrain: lambda u: sym(grad(u)),
mgis_bv.Hypothesis.Axisymmetrical: lambda r, u: sym(axi_grad(r, u))}
Grad = {mgis_bv.Hypothesis.Tridimensional: lambda u: grad(u),
mgis_bv.Hypothesis.PlaneStrain: lambda u: grad(u),
mgis_bv.Hypothesis.Axisymmetrical: lambda r, u: axi_grad(r, u)}
Id_Grad = {mgis_bv.Hypothesis.Tridimensional: lambda u: Identity(3)+grad(u),
mgis_bv.Hypothesis.PlaneStrain: lambda u: Identity(2)+grad(u),
mgis_bv.Hypothesis.Axisymmetrical: lambda r, u: Identity(3)+axi_grad(r, u)}
predefined_gradients = {
"Strain": SymGrad,
"TemperatureGradient": Grad,
"DeformationGradient": Id_Grad}
class MFrontNonlinearProblem(NonlinearProblem): class MFrontNonlinearProblem(NonlinearProblem):
def __init__(self, u, material, quadrature_degree=2): def __init__(self, u, material, quadrature_degree=2):
NonlinearProblem.__init__(self) NonlinearProblem.__init__(self)
...@@ -28,7 +43,7 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -28,7 +43,7 @@ class MFrontNonlinearProblem(NonlinearProblem):
self.ngauss = self.W0.dim() self.ngauss = self.W0.dim()
# Set data manager # Set data manager
self.material.set_data_manager(self.ngauss) self.material.set_data_manager(self.ngauss)
self.finite_strain = self.material.behaviour.getBehaviourType()=="StandardFiniteStrainBehaviour"
if self.material.hypothesis == mgis_bv.Hypothesis.Tridimensional: if self.material.hypothesis == mgis_bv.Hypothesis.Tridimensional:
assert self.u.geometric_dimension()==3, "Conflicting geometric dimension and material hypothesis" assert self.u.geometric_dimension()==3, "Conflicting geometric dimension and material hypothesis"
else: else:
...@@ -39,6 +54,7 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -39,6 +54,7 @@ class MFrontNonlinearProblem(NonlinearProblem):
"quadrature_scheme": "default"}) "quadrature_scheme": "default"})
if self.axisymmetric: if self.axisymmetric:
x = SpatialCoordinate(self.mesh) x = SpatialCoordinate(self.mesh)
self._r = x[0]
self.axi_coeff = 2*pi*abs(x[0]) self.axi_coeff = 2*pi*abs(x[0])
else: else:
self.axi_coeff = 1 self.axi_coeff = 1
...@@ -48,10 +64,25 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -48,10 +64,25 @@ class MFrontNonlinearProblem(NonlinearProblem):
self._Fext = None self._Fext = None
self._init = True self._init = True
self.state_variables = [] self.state_variables = dict.fromkeys(self.material.get_state_variable_names(), None)
self.gradients = dict.fromkeys(self.material.get_gradient_names(), None) self.gradients = dict.fromkeys(self.material.get_gradient_names(), None)
self.fluxes = {} self.fluxes = {}
self.predefined_registration()
def predefined_registration(self):
for g in self.gradients.keys():
if self.axisymmetric:
var = (self._r, self.u)
else:
var = (self.u,)
case = predefined_gradients.get(g, None)
if case is not None:
expr = case[self.material.hypothesis](*var)
self.register_gradient(g, expr)
print("Automatic registration of '{}' as {}.\n".format(g, str(expr)))
def define_form(self): def define_form(self):
# residual form (internal forces) # residual form (internal forces)
...@@ -60,13 +91,24 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -60,13 +91,24 @@ class MFrontNonlinearProblem(NonlinearProblem):
if self._Fext is not None: if self._Fext is not None:
self.L -= self._Fext self.L -= self._Fext
# tangent bilinear form # tangent bilinear form
try:
self.a = sum([inner(g.variation(self.u_), t*dg.variation(self.du)) self.a = sum([inner(g.variation(self.u_), t*dg.variation(self.du))
for (f, g) in zip(self.fluxes.values(), self.gradients.values()) for (f, g) in zip(self.fluxes.values(), self.gradients.values())
for (t, dg) in zip(f.tangent_blocks.values(), f.gradients)])*self.axi_coeff*self.dx for (t, dg) in zip(f.tangent_blocks.values(), f.gradients)])*self.axi_coeff*self.dx
except:
def register_gradient(self, name, expression, symmetric=None): print([(ufl.shape(g.variation(self.u_)), ufl.shape(t*dg.variation(self.du)))
for (f, g) in zip(self.fluxes.values(), self.gradients.values())
for (t, dg) in zip(f.tangent_blocks.values(), f.gradients)])
def register_gradient(self, name, expression):
pos = self.material.get_gradient_names().index(name) pos = self.material.get_gradient_names().index(name)
size = self.material.get_gradient_sizes()[pos] size = self.material.get_gradient_sizes()[pos]
vtype = self.material.behaviour.gradients[pos].type
if vtype == mgis_bv.VariableType.Tensor:
symmetric = False
elif vtype == mgis_bv.VariableType.Stensor:
symmetric = True
else:
symmetric = None
self.gradients.update({name: Gradient(self.u, expression, name, symmetric=symmetric)}) self.gradients.update({name: Gradient(self.u, expression, name, symmetric=symmetric)})
def set_loading(self, Fext): def set_loading(self, Fext):
...@@ -144,15 +186,19 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -144,15 +186,19 @@ class MFrontNonlinearProblem(NonlinearProblem):
mgis_bv.integrate(self.material.data_manager, mgis_bv.integrate(self.material.data_manager,
self.integration_type, 0, 0, self.material.data_manager.n) self.integration_type, 0, 0, self.material.data_manager.n)
if self.finite_strain: # if self.material.finite_strain:
local_project(self.strain_measure(self.u), self.Wsig, self.dx, self.strain) self.affect_gradients()
# copy the strain values to `MGIS` # else:
self.material.data_manager.s0.gradients[:, :] = self.strain.vector().get_local().reshape((self.material.data_manager.n, self.strain_dim))
else:
self.block_names = self.material.get_tangent_block_names() self.block_names = self.material.get_tangent_block_names()
self.block_shapes = self.material.get_tangent_block_sizes() self.block_shapes = self.material.get_tangent_block_sizes()
self.flattened_block_shapes = [s[0]*s[1] if len(s)==2 else s[0] for s in self.block_shapes] self.flattened_block_shapes = [s[0]*s[1] if len(s)==2 else s[0] for s in self.block_shapes]
self.update_tangent_blocks() self.update_tangent_blocks()
buff = 0
for (i, s) in enumerate(self.material.get_state_variable_names()):
state_var = self.state_variables[s]
block_shape = self.material.get_state_variable_sizes()[i]
self.material.data_manager.s0.internal_state_variables[:,buff:buff+block_shape] = state_var.vector().get_local().reshape(self.material.data_manager.n, block_shape)
buff += block_shape
def update_tangent_blocks(self): def update_tangent_blocks(self):
buff = 0 buff = 0
...@@ -160,6 +206,12 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -160,6 +206,12 @@ class MFrontNonlinearProblem(NonlinearProblem):
f, g = block f, g = block
t = self.fluxes[f].tangent_blocks[g] t = self.fluxes[f].tangent_blocks[g]
block_shape = self.flattened_block_shapes[i] block_shape = self.flattened_block_shapes[i]
if self.material.finite_strain and f == "Stress":
Ctv = t.vector().get_local()
mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager,
mgis_bv.FiniteStrainTangentOperator.DPK1_DF)
t.vector().set_local(Ctv)
else:
t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten()) t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten())
buff += block_shape buff += block_shape
...@@ -168,6 +220,12 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -168,6 +220,12 @@ class MFrontNonlinearProblem(NonlinearProblem):
for (i, f) in enumerate(self.material.get_flux_names()): for (i, f) in enumerate(self.material.get_flux_names()):
flux = self.fluxes[f] flux = self.fluxes[f]
block_shape = self.material.get_flux_sizes()[i] block_shape = self.material.get_flux_sizes()[i]
if self.material.finite_strain and f == "Stress":
pk1v = flux.function.vector().get_local()
mgis_bv.convertFiniteStrainStress(pk1v, self.material.data_manager,
mgis_bv.FiniteStrainStress.PK1)
flux.function.vector().set_local(pk1v)
else:
flux.function.vector().set_local(self.material.data_manager.s1.thermodynamic_forces[:,buff:buff+block_shape].flatten()) flux.function.vector().set_local(self.material.data_manager.s1.thermodynamic_forces[:,buff:buff+block_shape].flatten())
buff += block_shape buff += block_shape
...@@ -177,8 +235,34 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -177,8 +235,34 @@ class MFrontNonlinearProblem(NonlinearProblem):
gradient = self.gradients[g] gradient = self.gradients[g]
local_project(gradient(self.u), gradient.function_space, self.dx, gradient.function) local_project(gradient(self.u), gradient.function_space, self.dx, gradient.function)
block_shape = self.material.get_gradient_sizes()[i] block_shape = self.material.get_gradient_sizes()[i]
self.material.data_manager.s1.gradients[:, buff:buff+block_shape] = \ grad_vals = gradient.function.vector().get_local()
gradient.function.vector().get_local().reshape((self.material.data_manager.n, gradient.shape)) if gradient.shape > 0:
grad_vals = grad_vals.reshape((self.material.data_manager.n, gradient.shape))
else:
grad_vals = grad_vals[:, np.newaxis]
self.material.data_manager.s1.gradients[:, buff:buff+block_shape] = grad_vals
buff += block_shape
def affect_gradients(self):
buff = 0
for (i, g) in enumerate(self.material.get_gradient_names()):
gradient = self.gradients[g]
local_project(gradient(self.u), gradient.function_space, self.dx, gradient.function)
block_shape = self.material.get_gradient_sizes()[i]
grad_vals = gradient.function.vector().get_local()
if gradient.shape > 0:
grad_vals = grad_vals.reshape((self.material.data_manager.n, gradient.shape))
else:
grad_vals = grad_vals[:, np.newaxis]
self.material.data_manager.s0.gradients[:, buff:buff+block_shape] = grad_vals
buff += block_shape
def update_state_variables(self):
buff = 0
for (i, s) in enumerate(self.material.get_state_variable_names()):
state_var = self.state_variables[s]
block_shape = self.material.get_state_variable_sizes()[i]
state_var.vector().set_local(self.material.data_manager.s1.internal_state_variables[:,buff:buff+block_shape].flatten())
buff += block_shape buff += block_shape
...@@ -189,48 +273,40 @@ class MFrontNonlinearProblem(NonlinearProblem): ...@@ -189,48 +273,40 @@ class MFrontNonlinearProblem(NonlinearProblem):
0, 0, self.material.data_manager.n); 0, 0, self.material.data_manager.n);
# getting the stress and consistent tangent operator back to # getting the stress and consistent tangent operator back to
# the FEniCS world. # the FEniCS world.
if self.finite_strain: # if self.material.finite_strain:
pk1v = self.stress.vector().get_local() # pk1v = self.stress.vector().get_local()
mgis_bv.convertFiniteStrainStress(pk1v, self.material.data_manager, # mgis_bv.convertFiniteStrainStress(pk1v, self.material.data_manager,
mgis_bv.FiniteStrainStress.PK1) # mgis_bv.FiniteStrainStress.PK1)
self.stress.vector().set_local(pk1v) # self.stress.vector().set_local(pk1v)
Ctv = self.Ct.vector().get_local() # else:
mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager,
mgis_bv.FiniteStrainTangentOperator.DPK1_DF)
self.Ct.vector().set_local(Ctv)
else:
self.update_fluxes() self.update_fluxes()
self.update_tangent_blocks() self.update_tangent_blocks()
self.update_state_variables()
# sizes = self.material.get_state_variable_sizes()
# for i in range(len(self.state_variables)):
# size = sizes[i]
# s_values = self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten()
# state = self.state_variables[i]
# if isinstance(state, Flux):
# state.function.vector().set_local(s_values)
# shapes = [ufl.shape(t) for t in state.tangents]
# flattened_shapes = [s[0]*s[1] if len(s)==2 else s[0] if len(s)==1 else 1 for s in shapes]
# # FIXME: Assumes gradients of state variables are defined after fluxes and in right order
# for (j,t) in enumerate(state.tangents):
# s = flattened_shapes[j]
# t.vector().set_local(K[:,buff:buff+s].flatten())
# buff += s
# else:
# state.vector().set_local(s_values)
sizes = self.material.get_state_variable_sizes() def get_state_variable(self, name):
for i in range(len(self.state_variables)): if self.state_variables[name] is None:
size = sizes[i]
s_values = self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten()