diff --git a/src/materials/AbaqusMaterial.cpp b/src/materials/AbaqusMaterial.cpp new file mode 100755 index 0000000..35f63ad --- /dev/null +++ b/src/materials/AbaqusMaterial.cpp @@ -0,0 +1,329 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "util/Constants.h" +#include "AbaqusMaterial.h" + +using namespace jem; +using namespace jem::io; + +using jem::util::Properties; +using jive::Vector; + + +//======================================================================= +// class AbaqusMaterial +//======================================================================= + +//----------------------------------------------------------------------- +// static data +//----------------------------------------------------------------------- + +const char* AbaqusMaterial::UMAT_MATPROPS = "matProps"; +const char* AbaqusMaterial::UMAT_NSTATES = "nStates"; +const char* AbaqusMaterial::STORE_OLD_STRAIN = "storeOldStrain"; + +//----------------------------------------------------------------------- +// constructor & destructor +//----------------------------------------------------------------------- + + +AbaqusMaterial::AbaqusMaterial + + ( int rank, const Properties& globdat ) + : Material ( rank, globdat ), nStates_ ( 0 ) +{ + // Attempt to load the shared library 'libumat.so' + // RLTD_NOW resolves all symbols upon loading. + void* handle = dlopen(UMAT_SHARED_LIB, RTLD_NOW); + + // Check if 'handle' is a valid pointer + if ( handle == nullptr ) + { + std::cerr << "Error loading shared library: " + << dlerror() + << std::endl; + } + + // Function pointer to UMAT + umat = reinterpret_cast + ( dlsym( handle, "umat_" ) ); + + // Check function pointer validity + if ( umat == nullptr ) + { + std::cerr << "Error locating UMAT function: " + << dlerror() + << std::endl; + } + + // Abaqus strain vector is ordered differently from + // that of falcon. Therefore, perm_ is required to + // store the array permutations. + // Abaqus: [xx, yy, zz, xy, xz, yz] + // falcon: [xx, yy, zz, xy, yz, xz] + perm_.resize(6); + perm_ = {0,1,2,3,5,4}; + + // Set history allocation flag to false + allocd_ = false; +} + + +AbaqusMaterial::~AbaqusMaterial () +{} + + +//----------------------------------------------------------------------- +// configure +//----------------------------------------------------------------------- + + +void AbaqusMaterial::configure ( const Properties& props, + const Properties& globdat ) +{ + props.get ( matProps_, UMAT_MATPROPS ); + props.get ( nStates_, UMAT_NSTATES ); + props.get ( storeOld_, STORE_OLD_STRAIN ); +} + + +//----------------------------------------------------------------------- +// getConfig +//----------------------------------------------------------------------- + + +void AbaqusMaterial::getConfig ( const Properties& conf, + const Properties& globdat ) const +{ + conf.set ( UMAT_MATPROPS, matProps_ ); + conf.set ( UMAT_NSTATES, nStates_ ); + conf.set ( STORE_OLD_STRAIN, storeOld_ ); +} + + +//----------------------------------------------------------------------- +// updateConfig +//----------------------------------------------------------------------- + +void AbaqusMaterial::updateConfig () + +{ + // do nothing! +} + +//----------------------------------------------------------------------- +// update +//----------------------------------------------------------------------- + +void AbaqusMaterial::update + + ( Vector& stress, + Matrix& stiff, + const Vector& strain, + idx_t ipoint ) +{ + Vector dstrain; dstrain.resize(strain.size()); dstrain = 0.0; + + if ( allocd_ ) + { + dstrain = strain - preStrain_[ipoint].v; + } + + update ( stress, stiff, strain, dstrain, ipoint ); +} + +// ---------------------------------------------------------------- +// update (overloaded version) +// ---------------------------------------------------------------- + +void AbaqusMaterial::update + + ( Vector& stress, + Matrix& stiff, + const Vector& strain, + const Vector& dstrain, + int ipoint ) +{ + // Some constants + + const int strCount = STRAIN_COUNTS[rank_]; + const int matCount = matProps_.size(); + + // Initialize UMAT variables + + double STRESS[strCount] = {0.0}; // Zero vector + double STATEV[nStates_] = {0.0}; // Zero vector + double DDSDDE[strCount][strCount] = {0.0}; // Zero matrix + double SSE; // Not used + double SPD; // Not used + double SCD; // Not used + double RPL; // Not used + double DDSDDT; // Not used + double DRPLDE; // Not used + double DRPLDT; // Not used + double STRAN[strCount] = {0.0}; // Zero vector + double DSTRAN[strCount] = {0.0}; // Zero vector + double TIME; // Not used + double DTIME; // Not used + double TEMP; // Not used + double DTEMP; // Not used + double PREDEF; // Not used + double DPRED; // Not used + char CMNAME; // Not used + int NDI = 3; // 3 Direct stress + int NSHR = strCount - NDI; // Shear stress + int NTENS = strCount; // (= NDI + NSHR) + int NSTATV; // Not used int + double PROPS[matCount] = {0.0}; // Initialize mat props to zero + int NPROPS = matCount; // # material properties + double COORDS; // Not used + double DROT; // Not used + double PNEWDT; // Not used + double CELENT; // Not used + double DFGRD0; // Not used + double DFGRD1; // Not used + int NOEL; // Not used + int NPT; // Not used + int LAYER; // Not used + int KSPT; // Not used + int KSTEP; // Not used + int KIN; // Not used + + // Fill material properties + for ( int i = 0; i < matCount; i++ ) + { + PROPS[i] = matProps_[i]; + } + + // Fill strain and strain increment + for ( int i = 0; i < strCount; i++ ) + { + STRAN[i] = strain [perm_[i]]; + DSTRAN[i] = dstrain[perm_[i]]; + } + + // Fill state variables + if ( allocd_ ) + { + for ( int i = 0; i < nStates_; i++ ) + { + STATEV[i] = preState_[ipoint].v[i]; + } + } + + // Call Abaqus UMAT function + umat(&STRESS[0], &STATEV[0], &DDSDDE[0][0], &SSE, &SPD, + &SCD, &RPL, &DDSDDT, &DRPLDE, &DRPLDT, &STRAN[0], + &DSTRAN[0], &TIME, &DTIME, &TEMP, &DTEMP, &PREDEF, + &DPRED, &CMNAME, &NDI, &NSHR, &NTENS, &NSTATV, + &PROPS[0], &NPROPS, &COORDS, &DROT, &PNEWDT, &CELENT, + &DFGRD0, &DFGRD1, &NOEL, &NPT, &LAYER, &KSPT, &KSTEP, + &KIN); + + // Transfer Abaqus output to Jive variables (stress, stiff) + for ( int i = 0; i < strCount; i++ ) + { + stress[i] = STRESS[perm_[i]]; + + for ( int j = 0; j < strCount; j++ ) + { + stiff(i,j) = DDSDDE[perm_[i]][perm_[j]]; + } + } + + // Update state variables + if ( allocd_ ) + { + for ( int i = 0; i < nStates_; i++ ) + { + newState_[ipoint].v[i] = STATEV[i]; + } + } +} + + +//----------------------------------------------------------------------- +// allocPoints +//----------------------------------------------------------------------- + +void AbaqusMaterial::allocPoints ( int count ) +{ + preState_.reserve ( count ); + newState_.reserve ( count ); + + for ( idx_t ip = 0; ip < count; ++ip ) + { + preState_.pushBack ( Hist_( nStates_ ) ); + newState_.pushBack ( Hist_( nStates_ ) ); + } + + latestState_ = &preState_; + + if ( storeOld_ ) + { + preStrain_.reserve ( count ); + newStrain_.reserve ( count ); + + for ( idx_t ip = 0; ip < count; ++ip ) + { + preStrain_.pushBack ( Hist_( STRAIN_COUNTS[rank_] ) ); + newStrain_.pushBack ( Hist_( STRAIN_COUNTS[rank_] ) ); + } + + latestStrain_ = &preStrain_; + } + + // Set history allocation flag to true + allocd_ = true; +} + + +// -------------------------------------------------------------------- +// commit +// -------------------------------------------------------------------- + +void AbaqusMaterial::commit() +{ + newState_.swap ( preState_ ); + latestState_ = &( preState_ ); + + if ( storeOld_ ) + { + newStrain_.swap ( preStrain_ ); + latestStrain_ = &( preStrain_ ); + } +} + + +//----------------------------------------------------------------------- +// clone +//----------------------------------------------------------------------- + +Ref AbaqusMaterial::clone ( ) const + +{ + return newInstance ( *this ); +} + + +//----------------------------------------------------------------------- +// Hist_ constructor +//----------------------------------------------------------------------- + +AbaqusMaterial::Hist_::Hist_ ( int n ) +{ + v.resize ( n ); v = 0.0; +} \ No newline at end of file diff --git a/src/materials/AbaqusMaterial.h b/src/materials/AbaqusMaterial.h new file mode 100755 index 0000000..1fece63 --- /dev/null +++ b/src/materials/AbaqusMaterial.h @@ -0,0 +1,195 @@ +/** @file AbaqusMaterial.h + * @brief Abaqus UMAT interface. + * + * Author: R. Bharali, ritukesh.bharali@chalmers.se + * Date: 28 May 2024 + * + * Updates (when, what and who) + * - [XX YYYY 2024] + * + * @warning Not thoroughly tested. Not all UMAT + * function arguments are used. + */ + +#ifndef ABAQUSMATERIAL_H +#define ABAQUSMATERIAL_H + +#include +#include +#include +#include + +#include "Material.h" + +using jem::String; +using jem::util::Flex; +using jive::Vector; +using jive::IntVector; +using std::vector; + +// Define the Abaqus UMAT function signature +typedef void (*umatFunc)( + double* STRESS, // Stress + double* STATEV, // State variables + double* DDSDDE, // Material tangent + double* SSE, // Specific elastic energy + double* SPD, // Specific plastic dissipation + double* SCD, // Specific creep dissipation + double* RPL, // Volumetric heat generation per unit time + double* DDSDDT, // Tangent (stress w.r.t temperature) + double* DRPLDE, // Tangent (RPL w.r.t. strain) + double* DRPLDT, // Tangent (RPL w.r.t. temperature) + double* STRAN, // Strain + double* DSTRAN, // Strain increment + double* TIME, // Time + double* DTIME, // Time increment + double* TEMP, // Temperature + double* DTEMP, // Temperature increment + double* PREDEF, // Pre-defined field variables + double* DPRED, // Pre-defined field variables increment + char* CMNAME, // Material name + int* NDI, // # of direct stress components + int* NSHR, // # of shear stress components + int* NTENS, // Size of stress/strain array + int* NSTATV, // # of state variables + double* PROPS, // Material properties + int* NPROPS, // # of material properties + double* COORDS, // Coordinates of this point + double* DROT, // (3,3) rotation increment matrix + double* PNEWDT, // Ratio of new to old step-size + double* CELENT, // Characteristic element length + double* DFGRD0, // (3,3) Deformation grad at start of increment + double* DFGRD1, // (3,3) Deformation grad at end of increment + int* NOEL, // Element number + int* NPT, // Integration point number + int* LAYER, // Layer number + int* KSPT, // Section point number within the current layer + int* KSTEP, // Step number + int* KINC // Increment number + ); + +// Define the shared library as a macro +#define UMAT_SHARED_LIB "libumat.so" + +// ======================================================================= +// class AbaqusMaterial +// ======================================================================= + +/** @brief + * The class \c AbaqusMaterial implements an interface to Abaqus UMATs. + * In order to use this, first, the UMAT needs to be compiled as a + * shared library (libumat.so). Check out this Github repo + * to see how it is done. + * + * Below is an usage example, demonstrating how AbaqusMaterial is defined + * in the input file: + * + * \code + * material = + { + type = "AbaqusUMAT"; + matProps = [100., 0.2]; + nStates = 0; + storeOldStrain = true; + state = "PLANE_STRAIN"; // not used + rank = 2; // not used + }; + * \endcode + * + * Note that libumat.so must be presented in the same directory + * as the input file *.pro along with other input files. + * + * @warning Not thoroughly tested. Not all UMAT + * function arguments are used. + */ + + +class AbaqusMaterial : public Material +{ + public: + + typedef AbaqusMaterial Self; + + static const char* UMAT_MATPROPS; + static const char* UMAT_NSTATES; + static const char* STORE_OLD_STRAIN; + + explicit AbaqusMaterial + + ( int rank, + const Properties& globdat ); + + virtual void configure + + ( const Properties& props, + const Properties& globdat ); + + virtual void getConfig + + ( const Properties& conf, + const Properties& globdat ) const; + + virtual void updateConfig (); + + virtual void update + + ( Vector& stress, + Matrix& stiff, + const Vector& strain, + const idx_t ipoint ); + + virtual void update + + ( Vector& stress, + Matrix& stiff, + const Vector& strain, + const Vector& dstrain, + int ipoint ); + + virtual void allocPoints + + ( int ipCount ); + + virtual Ref clone () const; + + virtual void commit (); + + protected: + + virtual ~AbaqusMaterial (); + + private: + + Vector matProps_; //!< Abaqus UMAT material props + int nStates_; //!< Abaqus UMAT number of states + bool storeOld_; //!< Boolean flag to store old strain + umatFunc umat; //!< Function pointer to Abaqus UMAT + IntVector perm_; //!< Abaqus to Jive stress/strain permutations + + // Store current and old step states (internal variables) + + class Hist_ + { + public: + + Hist_( int n ); + + Vector v; + }; + + // Initialize state variables + + Flex preState_; + Flex newState_; + Flex* latestState_; + bool allocd_; + + // Initialize strains + + Flex preStrain_; + Flex newStrain_; + Flex* latestStrain_; + +}; + +#endif \ No newline at end of file diff --git a/src/materials/Material.cpp b/src/materials/Material.cpp index 543c037..25ee7ac 100644 --- a/src/materials/Material.cpp +++ b/src/materials/Material.cpp @@ -2,6 +2,7 @@ #include #include "Material.h" +#include "AbaqusMaterial.h" #include "HookeMaterial.h" #include "DamageMaterial.h" #include "AmorPhaseMaterial.h" @@ -240,6 +241,8 @@ Ref newMaterial mat = newInstance ( rank, globdat ); else if ( type == "MiehePhase" ) mat = newInstance ( rank, globdat ); + else if ( type == "AbaqusUMAT" ) + mat = newInstance ( rank, globdat ); else matProps.propertyError ( name, "Invalid material: " + type );