cl_Material_UserDefined.hpp

Go to the documentation of this file.

/*
 * BELFEM -- The Berkeley Lab Finite Element Framework
 * Copyright (c) 2026, The Regents of the University of California, through
 * Lawrence Berkeley National Laboratory (subject to receipt of any required
 * approvals from the U.S. Dept. of Energy).  All rights reserved.
 *
 * Developers: Christian Messe, Gregory Giard
 * 
 * See the top-level LICENSE file for the complete license and disclaimer.
 */
 
#ifndef BELFEM_CL_MATERIAL_USERDEFINED_HPP
#define BELFEM_CL_MATERIAL_USERDEFINED_HPP
 
#include "cl_Material.hpp"
 
// note: NEVER include cl_Vector.hpp or cl_Matrix.hpp here, or any class that uses it.
//       Doing so would break the API for the user defined materials.
 
namespace belfem
{
    namespace material
    {
//------------------------------------------------------------------------------
 
        template < typename T >
        T
        polyval( const std::vector< T > & coeffs, const T x )
        {
            const index_t n = coeffs.size() ;
            T y = coeffs[ 0 ];
 
            for( index_t k=1; k<n; ++k )
            {
                y *= x;
                y += coeffs[ k ];
            }
 
            return y;
        }
 
//------------------------------------------------------------------------------
 
        template < typename T >
        T
        dpolyval( const std::vector< T > & coeffs, const T x )
        {
            const index_t n = coeffs.size() - 1 ;
            T p = ( T ) n ;
            T dydx = p * coeffs[ 0 ];
 
            for( index_t k=1; k<n; ++k )
            {
                dydx *= x;
                p -= 1.0 ;
                dydx += p * coeffs[ k ];
            }
 
            return dydx;
        }
 
//------------------------------------------------------------------------------

        class UserDefinedMaterial : public Material
        {
            void *  mHandle ;  
 
            Cell< MatFunc1 * > mUserFunctions ;  
 
            Cell< Cell< real > > mPolynomials ;  
 
            real ( UserDefinedMaterial::*mMuFunction ) ( const real H, const real T ) const = nullptr ;  
 
            MatFunc2 * mUserMuFunction = nullptr ;  
            MatFunc3 * mUserRhoFunction = nullptr ;  
            MatFunc3 * mUserLambdaFunction = nullptr ;  
 
 
            public :

                UserDefinedMaterial( const string & aLibraryPath, const string & aLabel ) ;

                ~UserDefinedMaterial() override;

            void
            set_user_defined_function(
                const MaterialProperty   Property,
                const MaterialDependency Dependency,
                      MatFunc1 * Function ) override ;

            void
            set_user_defined_function(
                const MaterialProperty   Property,
                const MaterialDependency Dependency1,
                const MaterialDependency Dependency2,
                      MatFunc2 * Function ) override ;

            void
            set_user_defined_function(
                const MaterialProperty   Property,
                const MaterialDependency Dependency1,
                const MaterialDependency Dependency2,
                const MaterialDependency Dependency3,
                      MatFunc3 * Function ) override ;

            void
            set_user_defined_polynomial(
                const MaterialProperty Property,
                const Cell< real > & Coefficients ) override;
 
            void
            set_user_defined_polynomial(
               const MaterialProperty Property,
               const std::vector< real > & Coefficients ) override ;

            real
            evaluate_polynomial(
                const MaterialProperty  Property,
                const real T ) const override ;
 
            real
            evaluate_derivative_of_polynomial(
                const MaterialProperty  Property,
                const real T ) const override ;
 
        protected:

            real
            E_custom(const real T) const override
            {
                return mUserFunctions(  static_cast<size_t>(MaterialProperty::E) )( this, T );
            }

            real
            nu_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::nu))( this, T);
            }

            real
            cp_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::cp))( this, T);
            }

            real
            lambda_custom( const real T ) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::lambda))( this, T);
            }

            real
            lambda_custom(const real T, const real normB, const real angle) const override
            {
                return ( * mUserLambdaFunction ) (this, T, normB, angle );
            }

            real
            rho_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::rho))( this, T );
            }

            real
            rho_kohler(const real T, const real normB, const real angle) const override
            {
                return ( * mUserRhoFunction ) (this, T, normB, angle );
            }

            real
            jc_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::jc))( this, T );
            }

            real
            n_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::n))( this, T );
            }

            real
            alpha_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::alpha))( this, T);
            }

            real
            Rp02_custom(const real T) const override
            {
                return mUserFunctions( static_cast<size_t>(MaterialProperty::Rp02))( this, T);
            }

            real
            mu_custom( const real H, const real T ) const override
            {
                return ( this->*mMuFunction )( H, T );
            }
 
        private:

            real
            mu_polynomial( const real H, const real T ) const
            {
                return evaluate_polynomial( MaterialProperty::mu, T );
            }

            real
            mu_user( const real H, const real T ) const
            {
                return mUserMuFunction( this, H, T );
            }
 
            void
            set_user_defined_function( const MaterialProperty Property );
 
        };
 
//------------------------------------------------------------------------------
 
        inline real
        UserDefinedMaterial::evaluate_polynomial( const MaterialProperty  Property, const real T ) const
        {
            BELFEM_ASSERT( mPolynomials( static_cast<size_t>(Property) ).size() > 0,
                "return polynomial for %s of %s is not defined",
                to_string( Property ).c_str(),
                this->label().c_str()
            );
 
            return polyval( mPolynomials( static_cast<size_t>(Property) ).vector_data(), T );
        }
 
        inline real
        UserDefinedMaterial::evaluate_derivative_of_polynomial( const MaterialProperty  Property, const real T ) const
        {
            BELFEM_ASSERT( mPolynomials( static_cast<size_t>(Property) ).size() > 0,
                "return polynomial for %s of %s is not defined",
                to_string( Property ).c_str(),
                this->label().c_str()
            );
 
            return dpolyval( mPolynomials( static_cast<size_t>(Property) ).vector_data(), T );
        }
 
    }
}
#endif //BELFEM_CL_MATERIAL_USERDEFINED_HPP