Commit 65db9759 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Adds setup.py and operational wrapper

parent 88bca57f
A Python package for wrapping the usage of MFront in FEniCS
from dolfin import *
import mfront_wrapper as mf
import numpy as np
import ufl
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)
u = Function(V, name="Displacement")
def left(x, on_boundary):
return near(x[0], 0) and on_boundary
......@@ -26,17 +27,23 @@ 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 = XDMFFile("results/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")
material = mf.MFrontNonlinearMaterial('../materials/src/libBehaviour.so',
"LogarithmicStrainPlasticity")
problem = mf.MFrontNonlinearProblem(u, material)
problem.set_loading(dot(selfweight, u)*dx)
problem.bc = bc
p = problem.get_state_variable(name="EquivalentPlasticStrain")
assert (ufl.shape(p) == ())
epsel = problem.get_state_variable(name="ElasticStrain")
print(ufl.shape(epsel))
prm = problem.solver.parameters
prm["absolute_tolerance"] = 1e-6
prm["relative_tolerance"] = 1e-6
......
from dolfin import *
import mfront_wrapper as mf
import numpy as np
import ufl
hypothesis = "axisymmetric" # axisymmetric
......@@ -37,16 +38,18 @@ mat_prop = {"YoungModulus": E,
"PoissonRatio": nu,
"HardeningSlope": H,
"YieldStrength": sig0}
material = mf.NonlinearMaterial('materials/src/libBehaviour.so', 'IsotropicLinearHardeningPlasticity',
hypothesis=hypothesis,
material_properties=mat_prop)
material = mf.MFrontNonlinearMaterial('../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)
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")
file_results.parameters["flush_output"] = True
......
/*!
* \file TFEL/Material/IsotropicLinearHardeningPlasticity.hxx
* \brief this file implements the IsotropicLinearHardeningPlasticity Behaviour.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
......@@ -132,6 +132,7 @@ 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;
......@@ -448,23 +449,20 @@ template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const IsotropicLinearHardeningPlasticity<hypothesis,Type,false>& b)
{
using namespace std;
os << "eto : " << b.eto << '\n';
os << "deto : " << b.deto << '\n';
os << "sig : " << b.sig << '\n';
os << "dt : " << b.dt << endl;
os << "T : " << b.T << endl;
os << "dT : " << b.dT << endl;
os << "εᵗᵒ : " << b.eto << '\n';
os << "Δεᵗᵒ : " << b.deto << '\n';
os << "σ : " << b.sig << '\n';
os << "Δt : " << b.dt << '\n';
os << "young : " << b.young << '\n';
os << "nu : " << b.nu << '\n';
os << "H : " << b.H << '\n';
os << "s0 : " << b.s0 << '\n';
os << "eel : " << b.eel << '\n';
os << "deel : " << b.deel << '\n';
os << "Δeel : " << b.deel << '\n';
os << "p : " << b.p << '\n';
os << "dp : " << b.dp << '\n';
os << "Δp : " << b.dp << '\n';
os << "T : " << b.T << '\n';
os << "dT : " << b.dT << '\n';
os << "ΔT : " << b.dT << '\n';
os << "minimal_time_step_scaling_factor : " << b.minimal_time_step_scaling_factor << '\n';
os << "maximal_time_step_scaling_factor : " << b.maximal_time_step_scaling_factor << '\n';
return os;
......
/*!
* \file TFEL/Material/IsotropicLinearHardeningPlasticityBehaviourData.hxx
* \brief this file implements the IsotropicLinearHardeningPlasticityBehaviourData class.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
......@@ -85,6 +85,7 @@ 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;
......@@ -199,10 +200,8 @@ template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const IsotropicLinearHardeningPlasticityBehaviourData<hypothesis,Type,false>& b)
{
using namespace std;
os << "eto : " << b.eto << '\n';
os << "sig : " << b.sig << '\n';
os << "T : " << b.T << endl;
os << "εᵗᵒ : " << b.eto << '\n';
os << "σ : " << b.sig << '\n';
os << "young : " << b.young << '\n';
os << "nu : " << b.nu << '\n';
os << "H : " << b.H << '\n';
......
/*!
* \file TFEL/Material/IsotropicLinearHardeningPlasticityIntegrationData.hxx
* \brief this file implements the IsotropicLinearHardeningPlasticityIntegrationData class.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
......@@ -73,6 +73,7 @@ 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;
......@@ -168,11 +169,10 @@ template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const IsotropicLinearHardeningPlasticityIntegrationData<hypothesis,Type,false>& b)
{
using namespace std;
os << "deto : " << b.deto << '\n';
os << "dt : " << b.dt << endl;
os << "dT : " << b.dT << endl;
os << "dT : " << b.dT << '\n';
os << "Δεᵗᵒ : " << b.deto << '\n';
os << "σ : " << b.sig << '\n';
os << "Δt : " << b.dt << '\n';
os << "ΔT : " << b.dT << '\n';
return os;
}
......
/*!
* \file TFEL/Material/LogarithmicStrainPlasticity.hxx
* \brief this file implements the LogarithmicStrainPlasticity Behaviour.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
......@@ -142,6 +142,7 @@ 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;
......@@ -556,19 +557,16 @@ template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const LogarithmicStrainPlasticity<hypothesis,Type,false>& b)
{
using namespace std;
os << "eto : " << b.eto << '\n';
os << "deto : " << b.deto << '\n';
os << "sig : " << b.sig << '\n';
os << "dt : " << b.dt << endl;
os << "T : " << b.T << endl;
os << "dT : " << b.dT << endl;
os << "eel : " << b.eel << '\n';
os << "deel : " << b.deel << '\n';
os << "εᵗᵒ : " << b.eto << '\n';
os << "Δεᵗᵒ : " << b.deto << '\n';
os << "σ : " << b.sig << '\n';
os << "Δt : " << b.dt << '\n';
os << "εᵉˡ : " << b.eel << '\n';
os << "Δεᵉˡ : " << b.deel << '\n';
os << "p : " << b.p << '\n';
os << "dp : " << b.dp << '\n';
os << "Δp : " << b.dp << '\n';
os << "T : " << b.T << '\n';
os << "dT : " << b.dT << '\n';
os << "ΔT : " << b.dT << '\n';
os << "T_ : " << b.T_ << '\n';
os << "se : " << b.se << '\n';
os << "n : " << b.n << '\n';
......@@ -578,8 +576,8 @@ os << "s0 : " << b.s0 << '\n';
os << "Et : " << b.Et << '\n';
os << "minimal_time_step_scaling_factor : " << b.minimal_time_step_scaling_factor << '\n';
os << "maximal_time_step_scaling_factor : " << b.maximal_time_step_scaling_factor << '\n';
os << "theta : " << b.theta << '\n';
os << "epsilon : " << b.epsilon << '\n';
os << "θ : " << b.theta << '\n';
os << "ε : " << b.epsilon << '\n';
os << "iterMax : " << b.iterMax << '\n';
return os;
}
......
/*!
* \file TFEL/Material/LogarithmicStrainPlasticityBehaviourData.hxx
* \brief this file implements the LogarithmicStrainPlasticityBehaviourData class.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
......@@ -85,6 +85,7 @@ 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;
......@@ -177,11 +178,9 @@ template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const LogarithmicStrainPlasticityBehaviourData<hypothesis,Type,false>& b)
{
using namespace std;
os << "eto : " << b.eto << '\n';
os << "sig : " << b.sig << '\n';
os << "T : " << b.T << endl;
os << "eel : " << b.eel << '\n';
os << "εᵗᵒ : " << b.eto << '\n';
os << "σ : " << b.sig << '\n';
os << "εᵉˡ : " << b.eel << '\n';
os << "p : " << b.p << '\n';
os << "T : " << b.T << '\n';
return os;
......
/*!
* \file TFEL/Material/LogarithmicStrainPlasticityIntegrationData.hxx
* \brief this file implements the LogarithmicStrainPlasticityIntegrationData class.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
......@@ -73,6 +73,7 @@ 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;
......@@ -168,11 +169,10 @@ template<ModellingHypothesis::Hypothesis hypothesis,typename Type>
std::ostream&
operator <<(std::ostream& os,const LogarithmicStrainPlasticityIntegrationData<hypothesis,Type,false>& b)
{
using namespace std;
os << "deto : " << b.deto << '\n';
os << "dt : " << b.dt << endl;
os << "dT : " << b.dT << endl;
os << "dT : " << b.dT << '\n';
os << "Δεᵗᵒ : " << b.deto << '\n';
os << "σ : " << b.sig << '\n';
os << "Δt : " << b.dt << '\n';
os << "ΔT : " << b.dT << '\n';
return os;
}
......
......@@ -38,11 +38,14 @@ return policy;
extern "C"{
#endif /* __cplusplus */
MFRONT_SHAREDOBJ const char*
IsotropicLinearHardeningPlasticity_build_id = "";
MFRONT_SHAREDOBJ const char*
IsotropicLinearHardeningPlasticity_mfront_ept = "IsotropicLinearHardeningPlasticity";
MFRONT_SHAREDOBJ const char*
IsotropicLinearHardeningPlasticity_tfel_version = "3.3.0-dev";
IsotropicLinearHardeningPlasticity_tfel_version = "3.3.0";
MFRONT_SHAREDOBJ unsigned short IsotropicLinearHardeningPlasticity_mfront_mkt = 1u;
......@@ -106,7 +109,7 @@ MFRONT_SHAREDOBJ int IsotropicLinearHardeningPlasticity_ParametersTypes [] = {0,
MFRONT_SHAREDOBJ double IsotropicLinearHardeningPlasticity_minimal_time_step_scaling_factor_ParameterDefaultValue = 0.1;
MFRONT_SHAREDOBJ double IsotropicLinearHardeningPlasticity_maximal_time_step_scaling_factor_ParameterDefaultValue = 1.79769e+308;
MFRONT_SHAREDOBJ double IsotropicLinearHardeningPlasticity_maximal_time_step_scaling_factor_ParameterDefaultValue = 1.7976931348623e+308;
MFRONT_SHAREDOBJ unsigned short IsotropicLinearHardeningPlasticity_requiresStiffnessTensor = 0;
MFRONT_SHAREDOBJ unsigned short IsotropicLinearHardeningPlasticity_requiresThermalExpansionCoefficientTensor = 0;
......
/*!
* \file IsotropicLinearHardeningPlasticity.cxx
* \brief this file implements the IsotropicLinearHardeningPlasticity Behaviour.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Thomas Helfer
* \date 14 / 10 / 2016
*/
......
......@@ -43,11 +43,14 @@ return policy;
extern "C"{
#endif /* __cplusplus */
MFRONT_SHAREDOBJ const char*
LogarithmicStrainPlasticity_build_id = "";
MFRONT_SHAREDOBJ const char*
LogarithmicStrainPlasticity_mfront_ept = "LogarithmicStrainPlasticity";
MFRONT_SHAREDOBJ const char*
LogarithmicStrainPlasticity_tfel_version = "3.3.0-dev";
LogarithmicStrainPlasticity_tfel_version = "3.3.0";
MFRONT_SHAREDOBJ unsigned short LogarithmicStrainPlasticity_mfront_mkt = 1u;
......@@ -107,17 +110,17 @@ MFRONT_SHAREDOBJ const char * LogarithmicStrainPlasticity_Parameters[9] = {"Youn
"theta","epsilon","iterMax"};
MFRONT_SHAREDOBJ int LogarithmicStrainPlasticity_ParametersTypes [] = {0,0,0,0,0,0,0,0,2};
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_YoungModulus_ParameterDefaultValue = 2.1e+11;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_YoungModulus_ParameterDefaultValue = 210000000000;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_PoissonRatio_ParameterDefaultValue = 0;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_s0_ParameterDefaultValue = 2.5e+08;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_s0_ParameterDefaultValue = 250000000;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_Et_ParameterDefaultValue = 1e+06;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_Et_ParameterDefaultValue = 1000000;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_minimal_time_step_scaling_factor_ParameterDefaultValue = 0.1;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_maximal_time_step_scaling_factor_ParameterDefaultValue = 1.79769e+308;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_maximal_time_step_scaling_factor_ParameterDefaultValue = 1.7976931348623e+308;
MFRONT_SHAREDOBJ double LogarithmicStrainPlasticity_theta_ParameterDefaultValue = 1;
......
/*!
* \file LogarithmicStrainPlasticity.cxx
* \brief this file implements the LogarithmicStrainPlasticity Behaviour.
* File generated by tfel version 3.3.0-dev
* File generated by tfel version 3.3.0
* \author Helfer Thomas
* \date 5 / 12 / 13
*/
......
This diff is collapsed.
# Makefile generated by mfront.
# tfel
# Version : 3.3.0-dev
# Version : 3.3.0
# Compiled with on -
# Please submit bug at tfel-contact@cea.fr
......@@ -22,7 +22,7 @@ makefiles = $(makefiles2)
all : libBehaviour.so
libBehaviour.so : LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o
libBehaviour.so : IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o
@$(CXX) -shared $^ -o $@ -L"$(strip $(shell tfel-config --library-path))" $(patsubst %,-l%,$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants))
install :
......
......@@ -6,10 +6,10 @@ prefix : "lib";
suffix : "so";
install_path : "";
sources : {
"LogarithmicStrainPlasticity-generic.cxx",
"LogarithmicStrainPlasticity.cxx",
"IsotropicLinearHardeningPlasticity-generic.cxx",
"IsotropicLinearHardeningPlasticity.cxx"
"IsotropicLinearHardeningPlasticity.cxx",
"LogarithmicStrainPlasticity-generic.cxx",
"LogarithmicStrainPlasticity.cxx"
};
cppflags : {
"$(shell tfel-config --cppflags --compiler-flags)"
......@@ -24,26 +24,26 @@ link_libraries : {
"$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants)"
};
epts : {
"LogarithmicStrainPlasticity_AxisymmetricalGeneralisedPlaneStrain",
"LogarithmicStrainPlasticity_Axisymmetrical",
"LogarithmicStrainPlasticity_PlaneStrain",
"LogarithmicStrainPlasticity_GeneralisedPlaneStrain",
"LogarithmicStrainPlasticity_Tridimensional",
"IsotropicLinearHardeningPlasticity_AxisymmetricalGeneralisedPlaneStrain",
"IsotropicLinearHardeningPlasticity_Axisymmetrical",
"IsotropicLinearHardeningPlasticity_PlaneStrain",
"IsotropicLinearHardeningPlasticity_GeneralisedPlaneStrain",
"IsotropicLinearHardeningPlasticity_Tridimensional"
"IsotropicLinearHardeningPlasticity_Tridimensional",
"LogarithmicStrainPlasticity_AxisymmetricalGeneralisedPlaneStrain",
"LogarithmicStrainPlasticity_Axisymmetrical",
"LogarithmicStrainPlasticity_PlaneStrain",
"LogarithmicStrainPlasticity_GeneralisedPlaneStrain",
"LogarithmicStrainPlasticity_Tridimensional"
};
};
headers : {
"MFront/GenericBehaviour/LogarithmicStrainPlasticity-generic.hxx",
"TFEL/Material/LogarithmicStrainPlasticity.hxx",
"TFEL/Material/LogarithmicStrainPlasticityBehaviourData.hxx",
"TFEL/Material/LogarithmicStrainPlasticityIntegrationData.hxx",
"MFront/GenericBehaviour/IsotropicLinearHardeningPlasticity-generic.hxx",
"TFEL/Material/IsotropicLinearHardeningPlasticity.hxx",
"TFEL/Material/IsotropicLinearHardeningPlasticityBehaviourData.hxx",
"TFEL/Material/IsotropicLinearHardeningPlasticityIntegrationData.hxx"
"TFEL/Material/IsotropicLinearHardeningPlasticityIntegrationData.hxx",
"MFront/GenericBehaviour/LogarithmicStrainPlasticity-generic.hxx",
"TFEL/Material/LogarithmicStrainPlasticity.hxx",
"TFEL/Material/LogarithmicStrainPlasticityBehaviourData.hxx",
"TFEL/Material/LogarithmicStrainPlasticityIntegrationData.hxx"
};
};
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 21 16:33:46 2020
@author: bleyerj
"""
from dolfin import *
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
from .nonlinear_material import MFrontNonlinearMaterial
from .nonlinear_problem import MFrontNonlinearProblem
\ No newline at end of file
from dolfin import *
import ufl
import mgis.behaviour as mgis_bv
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
mgis_hypothesis = {"plane_strain": mgis_bv.Hypothesis.PlaneStrain,
"plane_stress": mgis_bv.Hypothesis.PlaneStress,
"3d": mgis_bv.Hypothesis.Tridimensional,
"axisymmetric": mgis_bv.Hypothesis.Axisymmetrical}
def local_project(v, V, dx, u=None):
"""
projects v on V with custom quadrature scheme dedicated to
FunctionSpaces V of `Quadrature` type
if u is provided, result is appended to u
"""
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dx
b_proj = inner(v, v_)*dx
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
def symmetric_tensor_to_vector(T, T22=0):
""" Return symmetric tensor components in vector form notation following MFront conventions
T22 can be specified when T is only (2,2)"""
if ufl.shape(T)==(2, 2):
return as_vector([T[0, 0], T[1, 1], T22, sqrt(2)*T[0, 1]])
elif ufl.shape(T)==(3, 3):
return as_vector([T[0, 0], T[1, 1], T[2, 2], sqrt(2)*T[0, 1], sqrt(2)*T[0, 2], sqrt(2)*T[1, 2]])
else:
raise NotImplementedError
def nonsymmetric_tensor_to_vector(T, T22=0):
""" Return nonsymmetric tensor components in vector form notation following MFront conventions
T22 can be specified when T is only (2,2) """
if ufl.shape(T)==(2, 2):
return as_vector([T[0, 0], T[1, 1], T22, T[0, 1], T[1, 0]])
elif ufl.shape(T)==(3, 3):
return as_vector([T[0, 0], T[1, 1], T[2, 2], T[0, 1], T[1, 0], T[0, 2], T[2, 0], T[1, 2], T[2, 1]])
else:
raise NotImplementedError
def axi_grad(r, v):
"""
Axisymmetric gradient in cylindrical coordinate (er, etheta, ez) for:
* a scalar v(r, z)
* a 2d-vectorial (vr(r,z), vz(r, z))
* a 3d-vectorial (vr(r,z), 0, vz(r, z))
"""
if ufl.shape(v)==(3,):
return as_matrix([[v[0].dx(0), -v[1]/r, v[0].dx(1)],
[v[1].dx(0), v[0]/r, v[1].dx(1)],
[v[2].dx(0), 0, v[2].dx(1)]])
elif ufl.shape(v)==(2,):
return as_matrix([[v[0].dx(0), v[0].dx(1), 0],
[v[1].dx(0), v[1].dx(1), 0],
[0, 0, v[0]/r]])
elif ufl.shape(v)==():
return as_vector([v.dx(0), 0, v.dx(1)])
else:
raise NotImplementedError
def symmetric_gradient(g):
""" Return symmetric gradient components in vector form"""
return symmetric_tensor_to_vector(sym(g))
def transformation_gradient(g):
""" Return transformation gradient components in vector form"""
return nonsymmetric_tensor_to_vector(Identity(v.geometric_dimension())+g, T22=1)
def gradient(g):
""" Return displacement gradient components in vector form"""
return nonsymmetric_tensor_to_vector(g)
def get_quadrature_element(cell, degree, dim=0):
if dim == 0 or dim == ():
return FiniteElement("Quadrature", cell, degree=degree,
quad_scheme='default')
elif type(dim) == int or len(dim)==1:
return VectorElement("Quadrature", cell, degree=degree,
dim=dim, quad_scheme='default')
elif type(dim) == tuple:
return TensorElement("Quadrature", cell, degree=degree, shape=dim, quad_scheme='default')
else:
raise ValueError("Wrong shape for dim=", dim)