Skip to content

Commit

Permalink
Add an optional temperature field to the umat wrapper
Browse files Browse the repository at this point in the history
  • Loading branch information
@pruliere
pruliere committed Jan 30, 2024
1 parent 2f04f56 commit ca36975
Show file tree
Hide file tree
Showing 3 changed files with 281 additions and 8 deletions.
4 changes: 3 additions & 1 deletion ...oon-python-builder/include/simcoon/python_wrappers/Libraries/Continuum_mechanics/umat.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
#pragma once
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
namespace py=pybind11;
using namespace std;

namespace simpy {
pybind11::tuple launch_umat(const std::string& umat_name_py, const pybind11::array_t<double> &etot_py, const pybind11::array_t<double> &Detot_py, const pybind11::array_t<double> &sigma_py, const pybind11::array_t<double> &DR_py, const pybind11::array_t<double> &props_py, const pybind11::array_t<double> &statev_py, const float Time, const float DTime, const pybind11::array_t<double> &Wm_py);
py::tuple launch_umat(const std::string& umat_name_py, const py::array_t<double> &etot_py, const py::array_t<double> &Detot_py, const py::array_t<double> &sigma_py, const py::array_t<double> &DR_py, const py::array_t<double> &props_py, const py::array_t<double> &statev_py, const float Time, const float DTime, const py::array_t<double> &Wm_py, const std::optional<py::array_t<double>> &T_py);
}
281 changes: 276 additions & 5 deletions simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/umat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ namespace py=pybind11;

namespace simpy {

py::tuple launch_umat(const std::string& umat_name_py, const py::array_t<double> &etot_py, const py::array_t<double> &Detot_py, const py::array_t<double> &sigma_py, const py::array_t<double> &DR_py, const py::array_t<double> &props_py, const py::array_t<double> &statev_py, const float Time, const float DTime, const py::array_t<double> &Wm_py){
py::tuple launch_umat(const std::string& umat_name_py, const py::array_t<double> &etot_py, const py::array_t<double> &Detot_py, const py::array_t<double> &sigma_py, const py::array_t<double> &DR_py, const py::array_t<double> &props_py, const py::array_t<double> &statev_py, const float Time, const float DTime, const py::array_t<double> &Wm_py, const std::optional<py::array_t<double>> &T_py){
//Get the id of umat
std::map<string, int> list_umat;
list_umat = { {"UMEXT",0},{"UMABA",1},{"ELISO",2},{"ELIST",3},{"ELORT",4},{"EPICP",5},{"EPKCP",6},{"EPCHA",7},{"SMAUT",8},{"LLDM0",9},{"ZENER",10},{"ZENNK",11},{"PRONK",12},{"EPHIC",17},{"EPHIN",18},{"SMAMO",19},{"SMAMC",20},{"MIHEN",100},{"MIMTN",101},{"MISCN",103},{"MIPLN",104} };
Expand Down Expand Up @@ -99,7 +99,6 @@ namespace simpy {
//simcoon::umat_plasticity_chaboche_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt);
break;
}

case 8: {
umat_function_2 = &simcoon::umat_sma_unified_T;
arguments_type = 2;
Expand Down Expand Up @@ -171,13 +170,15 @@ namespace simpy {
const int ncomp=ndi+nshr;
const bool start = false;
double tnew_dt = 0;//usefull ?
double T = 0; double DT = 0; //modify to let the program set the actual temperature
double T = 0; //default value
double DT = 0; //not used. T should be the current temperature (=T+DT)
//bool use_temp;
//if (T.n_elem == 0.) use_temp = false;
//else use_temp = true;

//py::print("ndim = ", etot_py.ndim());
if (etot_py.ndim() == 1) {
if (etot_py.ndim() == 1) {
if (T_py.has_value()) T=T_py.value().at(0);
vec etot = carma::arr_to_col_view(etot_py);
vec Detot = carma::arr_to_col_view(Detot_py);
vec sigma = carma::arr_to_col(sigma_py); //copy data because values are changed by the umat and returned to python
Expand Down Expand Up @@ -210,7 +211,15 @@ namespace simpy {
return py::make_tuple(carma::col_to_arr(sigma, false), carma::col_to_arr(statev, false), carma::col_to_arr(Wm, false), carma::mat_to_arr(Lt, false));
}
else if (etot_py.ndim() == 2) {

bool use_temp;
vec vec_T;
if (T_py.has_value()) {
vec_T = carma::arr_to_col_view(T_py.value());
use_temp = true;
}
else {
use_temp = false;
}
mat etot = carma::arr_to_mat_view(etot_py);
int nb_points = etot.n_cols; //number of material points
mat Detot = carma::arr_to_mat_view(Detot_py);
Expand All @@ -232,6 +241,7 @@ namespace simpy {
if (pt < list_props.n_cols) props = list_props.col(pt); //if list_props has only one element, we keep only this one (assuming homogeneous material)
vec statev = list_statev.unsafe_col(pt);
vec sigma = list_sigma.unsafe_col(pt);
if (use_temp) T=vec_T(pt);

switch (arguments_type) {

Expand All @@ -257,3 +267,264 @@ namespace simpy {



/* py::tuple launch_umat_T(const std::string& umat_name_py, const py::array_t<double> &etot_py, const py::array_t<double> &Detot_py, const py::array_t<double> &sigma_py, const py::array_t<double> &DR_py, const py::array_t<double> &props_py, const py::array_t<double> &statev_py, const float Time, const float DTime, const py::array_t<double> &Wm_py, py::array_t<double> &T){
//Get the id of umat
std::map<string, int> list_umat;
list_umat = { {"UMEXT",0},{"UMABA",1},{"ELISO",2},{"ELIST",3},{"ELORT",4},{"EPICP",5},{"EPKCP",6},{"EPCHA",7},{"SMAUT",8},{"LLDM0",9},{"ZENER",10},{"ZENNK",11},{"PRONK",12},{"EPHIC",17},{"EPHIN",18},{"SMAMO",19},{"SMAMC",20},{"MIHEN",100},{"MIMTN",101},{"MISCN",103},{"MIPLN",104} };
int id_umat = list_umat[umat_name_py];
int arguments_type; //depends on the argument used in the umat
void (*umat_function)(const arma::vec &, const arma::vec &, arma::vec &, double &, arma::mat &, arma::mat &, arma::mat &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &);
switch (id_umat) {
case 2: {
umat_function = &simcoon::umat_elasticity_iso_T;
arguments_type = 1;
//umat_elasticity_iso_T(umat_T->Etot, umat_T->DEtot, umat_T->sigma, umat_T->r, umat_T->dSdE, umat_T->dSdT, umat_T->drdE, umat_T->drdT, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_T->nstatev, umat_T->statev, umat_T->T, umat_T->DT, Time, DTime, umat_T->Wm(0), umat_T->Wm(1), umat_T->Wm(2), umat_T->Wm(3), umat_T->Wt(0), umat_T->Wt(1), umat_T->Wt(2), ndi, nshr, start, tnew_dt);
break;
}
case 3: {
umat_function = &simcoon::umat_elasticity_trans_iso_T;
arguments_type = 4;
break;
}
case 4: {
umat_function = &simcoon::umat_elasticity_ortho_T;
arguments_type = 4;
break;
}
/* case 5: {
umat_function = &simcoon::umat_plasticity_iso_CCP;
arguments_type = 1;
//simcoon::umat_plasticity_iso_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt);
break;
}
case 6: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function = &simcoon::umat_plasticity_kin_iso_CCP;
arguments_type = 1;
//simcoon::umat_plasticity_kin_iso_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt);
}
break;
}
case 7: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function = &simcoon::umat_plasticity_chaboche_CCP;
arguments_type = 1;
//simcoon::umat_plasticity_chaboche_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt);
}
break;
}
case 8: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_2 = &simcoon::umat_sma_unified_T;
arguments_type = 2;
//simcoon::umat_sma_unified_T(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 9: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function = &simcoon::umat_damage_LLD_0;
arguments_type = 1;
//simcoon::umat_damage_LLD_0(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt);
}
break;
}
case 10: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_2 = &simcoon::umat_zener_fast;
arguments_type = 2;
//simcoon::umat_zener_fast(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 11: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_2 = &simcoon::umat_zener_Nfast;
arguments_type = 2;
//simcoon::umat_zener_Nfast(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 12: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_2 = &simcoon::umat_prony_Nfast;
arguments_type = 2;
//simcoon::umat_prony_Nfast(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 17: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_2 = &simcoon::umat_plasticity_hill_isoh_CCP;
arguments_type = 2;
//simcoon::umat_plasticity_hill_isoh_CCP(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 18: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_2 = &simcoon::umat_plasticity_hill_isoh_CCP_N;
arguments_type = 2;
//simcoon::umat_plasticity_hill_isoh_CCP_N(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 19: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_3 = &simcoon::umat_sma_mono;
arguments_type = 3;
//simcoon::umat_sma_mono(etot, Detot, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
case 20: {
if (thermomechancial){
arguments_type = 0;
}
else {
umat_function_3 = &simcoon::umat_sma_mono_cubic;
arguments_type = 3;
//simcoon::umat_sma_mono_cubic(etot, Detot, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt);
}
break;
}
default: {
//py::print("Error: The choice of Umat could not be found in the umat library \n");
throw std::invalid_argument( "The choice of Umat could not be found in the umat library." );
//exit(0);
}
}
if (arguments_type == 0) throw std::invalid_argument( "The choice of Umat could not be found in the umat library. Check the thermomecanical argument." );
//py::print("id_umat ok");
//scalar needed to launch umat
const int solver_type = 0;
const int ndi=3;
const int nshr=3;
const int ncomp=ndi+nshr;
const bool start = false;
double tnew_dt = 0;//usefull ?
double T = 0; double DT = 0; //modify to let the program set the actual temperature
//bool use_temp;
//if (T.n_elem == 0.) use_temp = false;
//else use_temp = true;
//py::print("ndim = ", etot_py.ndim());
if (etot_py.ndim() == 1) {
vec etot = carma::arr_to_col_view(etot_py);
vec Detot = carma::arr_to_col_view(Detot_py);
vec sigma = carma::arr_to_col(sigma_py); //copy data because values are changed by the umat and returned to python
mat DR = carma::arr_to_mat_view(DR_py);
vec props = carma::arr_to_col_view(props_py);
vec statev = carma::arr_to_col(statev_py); //copy data because values are changed by the umat and returned to python
vec Wm = carma::arr_to_col(Wm_py); //copy data because values are changed by the umat and returned to python
mat L(ncomp, ncomp);
mat Lt(ncomp, ncomp);
vec sigma_in = zeros(1); //not used
int nprops = props.n_elem;
int nstatev = statev.n_elem;
switch (arguments_type) {
case 1: {
umat_function(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, solver_type, tnew_dt);
break;
}
case 2: {
umat_function_2(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt);
break;
}
case 3: {
umat_function_3(etot, Detot, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt);
break;
}
case 4: {
umat_elasticity_iso_T(etot, Detot, sigma, umat_T->r, umat_T->dSdE, umat_T->dSdT, umat_T->drdE, umat_T->drdT, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), umat_T->Wt(0), umat_T->Wt(1), umat_T->Wt(2), ndi, nshr, start, tnew_dt);
break;
}
}
//py::print("umat_done");
return py::make_tuple(carma::col_to_arr(sigma, false), carma::col_to_arr(statev, false), carma::col_to_arr(Wm, false), carma::mat_to_arr(Lt, false));
}
else if (etot_py.ndim() == 2) {
mat etot = carma::arr_to_mat_view(etot_py);
int nb_points = etot.n_cols; //number of material points
mat Detot = carma::arr_to_mat_view(Detot_py);
mat list_sigma = carma::arr_to_mat(sigma_py); //copy data because values are changed by the umat and returned to python
cube DR = carma::arr_to_cube_view(DR_py);
mat list_props = carma::arr_to_mat_view(props_py);
mat list_statev = carma::arr_to_mat(statev_py); //copy data because values are changed by the umat and returned to python
mat Wm = carma::arr_to_mat(Wm_py); //copy data because values are changed by the umat and returned to python
cube L(ncomp, ncomp, nb_points);
cube Lt(ncomp, ncomp, nb_points);
vec sigma_in = zeros(1); //not used
int nprops = list_props.n_rows;
int nstatev = list_statev.n_rows;
vec props(ncomp);
for (int pt = 0; pt < nb_points; pt++) {
//if (use_temp) T = list_T(pt);
if (pt < list_props.n_cols) props = list_props.col(pt); //if list_props has only one element, we keep only this one (assuming homogeneous material)
vec statev = list_statev.unsafe_col(pt);
vec sigma = list_sigma.unsafe_col(pt);
switch (arguments_type) {
case 1: {
umat_function(etot.col(pt), Detot.col(pt), sigma, Lt.slice(pt), L.slice(pt), sigma_in, DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0,pt), Wm(1,pt), Wm(2,pt), Wm(3,pt), ndi, nshr, start, solver_type, tnew_dt);
break;
}
case 2: {
umat_function_2(etot.col(pt), Detot.col(pt), sigma, Lt.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0,pt), Wm(1,pt), Wm(2,pt), Wm(3,pt), ndi, nshr, start, tnew_dt);
break;
}
case 3: {
umat_function_3(etot.col(pt), Detot.col(pt), sigma, Lt.slice(pt), L.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0,pt), Wm(1,pt), Wm(2,pt), Wm(3,pt), ndi, nshr, start, tnew_dt);
break;
}
}
}
return py::make_tuple(carma::mat_to_arr(list_sigma, false), carma::mat_to_arr(list_statev, false), carma::mat_to_arr(Wm, false), carma::cube_to_arr(Lt, false));
}
}
}
*/
4 changes: 2 additions & 2 deletions simcoon-python-builder/src/python_wrappers/boostpython_smartplus_wrappers.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

#include <CGAL/Simple_cartesian.h>

Expand Down Expand Up @@ -230,8 +231,7 @@ PYBIND11_MODULE(simmit, m) {
m.def("stress_convert", &stress_convert, "sigma"_a, "F"_a, "converter_key"_a, "J"_a=0., "copy"_a=true, "Provides the first Piola Kirchoff stress tensor from the Cauchy stress tensor");

//umat
m.def("umat", &launch_umat);
//py::tuple umat(const std::string& umat_name_py, const py::array_t<double> &etot_py, const py::array_t<double> &Detot_py, const py::array_t<double> &sigma_py, const py::array_t<double> &DR_py, const py::array_t<double> &props_py, const py::array_t<double> &statev_py, const float Time, const float DTime, const py::array_t<double> &Wm_py){
m.def("umat", &launch_umat, "umat_name"_a, "etot"_a, "Detot"_a, "sigma"_a, "DR"_a, "props"_a, "statev"_a, "time"_a, "dtime"_a, "Wm"_a, "temp"_a = py::none());

// Register the from-python converters for read and solver
m.def("read_matprops", &read_matprops);
Expand Down

0 comments on commit ca36975

Please sign in to comment.