Implementation

Utility Structures

A few structures/classes are provided to help in writing and communicating with constitutive models. Some are provided by “base007.dll”, which defines the base level of functionality common to all plugin interfaces. Others are provided by “conmodel007.dll”, which defines the specific interface used by the constitutive model system. In all cases, full documentation of the class definition is available in the base interface module of the Programmer’s Interface documentation.

“base007.dll” provides:

Int, UInt, Byte, UByte, Double, Float, etc. — defined in “base/src/basedef.h”. These types are substitutions for the standard C++ types int, unsigned, char, double, etc. These types are used instead to define consistent size definitions. For instance, an Int is guaranteed to be asigned 32-bit quantity regardless of the platform.

String — defined in “base/src/basebool.h”. This class is derived from std::basic_string<wchar_t>, or the ANSI C++ standard string class. The String class adds a few handy utility functions, such as string, to numeric conversion.

DVect3 — defined in “base/src/vect.h”. This class is actually the double instance of the template class Vector3<T>. Similar predefined types are IVect (Vector3<Int>) and UVect (Vector3<UInt>). This class allows one to treat a three-dimensional vector as a primitive type, with full operator overloading for convenient syntax.

Variant — defined in “base/src/variant.h”. This defines a type that can represent many different primitive types (e.g., String, Double, Int, etc.). This class is used to pass properties to and from the constitutive model.

Axes — defined in “base/src/axes.h”. This class allows for the definition of an orthonormal basis. This basis can be used to convert coordinates to and from a “global” basis represented by the traditional Cartesian axes.

SymTensor — defined in “base/src/symtensor.h”. This class defines a 3 \(\times\) 3 symmetric tensor, typically used to represent stresses and strains. Member access is available through the s11(), s12(), s13(), s22(), etc. functions. Member modification is available through the rs11(), rs12(), etc. functions. In addition, eigenvector information (or principal directions and values) can be obtained through the getEigenInfo() function. The helper class SymTensorInfo is used to allow the user to modify principal values while maintaining principal directions, and build up a new SymTensor with the result.

Orientation3 - defined in “base/src/orientation.h”. This class provides storage and manipulation of an orientation, or a direction in space. This orientation can be specified either by a normal vector, or by a dip and dip direction.

In addition to the ConstitutiveModel and State interfaces, the “conmodel007.dll” provides the following two utility functions in “convert.h”: getYPFromBS() and getBSfromYP(). These functions can be used to convert Young’s modulus and Poisson’s ratio values to bulk and shear modulus values, and vice versa.

Example Constitutive Model

The source codes of all constitutive models used by FLAC3D are provided for the user to inspect or adapt. Here we extract, for illustration, parts of the Mohr-Coulomb elastic/plastic model contained in files “modelexample.*”. modelexample.h provides the class specification for the model, which also includes the definition of the model’s unique type number. Note that there are more protected variables than property names in the ModelElastic base class (see the getProperties() member function). In this model, some of the variables are for internal use only: they occupy memory in each zone, but they are not available for the user of FLAC3D to change or print out. Also note that the getProperty()/setProperty interface is used for Save/Restore.

modelexample.h:

#pragma once

#include "models/elastic/src/modelelastic.h"

namespace models {
    class ModelExample : public ModelElastic {
    public:
        ModelExample(unsigned short option=0);
        virtual String        getName() const override;
        virtual String        getFullName() const override;
        virtual UInt          getMinorVersion() const override; 
        virtual String        getProperties() const override;
        virtual String        getStates() const override;
        virtual Variant       getProperty(UInt index) const override;
        virtual void          setProperty(UInt index,const Variant &p,UInt restoreVersion=0) override;
        virtual ModelExample *clone() const override { return new ModelExample(); }      
        virtual void          copy(const ConstitutiveModel *mod) override;
        virtual void          run(UByte dim,State *s) override; 
        virtual void          initialize(UByte dim,State *s) override;
        // Optional
        virtual bool          isPropertyAdvanced(UInt i) const override;
        virtual bool          supportsStressStrengthRatio() const override { return true; }
        virtual bool          supportsPropertyScaling() const override { return true; }
        virtual Double        getStressStrengthRatio(const SymTensor &st) const override;
        virtual void          scaleProperties(const Double &scale,const std::vector<UInt> &props) override;

    private:
        virtual bool          updateParameters() override;
        virtual bool          updateParameters(bool bEUpdated,Double *sf1=nullptr,Double *sf3=nullptr);
        virtual Double        moduliReduction(const Double &factor) override;
        virtual void          apexCorrection(const Double &fric,State *s,DVect3 *prin,UInt *iPlasticity=nullptr,bool bBrittle=false);
        virtual void          tensionCorrection(State *s,DVect3 *prin,UInt *iPlasticity,const Double &ftz,bool bBrittle=false);
        virtual void          shearCorrection(State *s,DVect3 *prin,UInt *iPlasticity,const Double &fs);

        Double cohesion_ = 0.0;
        Double friction_ = 0.0;
        Double dilation_ = 0.0;
        Double tension_ = 0.0;
        bool brittle_ = false;
        Double nph_  = 0.0;
        Double csn_ = 0.0;
        Double nps_ = 0.0;
        Double rc_ = 0.0;
        Double sc1_ = 0.0;
        Double sc2_ = 0.0;
        Double sc3_ = 0.0;
    };
} // namespace models

// EOF

The constructor for this model was listed in Typical Model Constructor. The example initialize and run functions provides listings of the member functions for initialization and execution (“running’‘). Note that, to save time, protected model variables e1_, e2_, g2_, etc. are not computed at each cycle, rather during initializatin in the ModelElastic::updateParameters method. Also note the use of the State structure in providing strain increments and stresses. Separate sections could be provided in every model for execution in two and three dimensions, to allow the same models to be used efficiently in FLAC or UDEC. In this example, the 2D section is identical to the 3D section. Please refer to the file “modelexample.cpp” for listings of member functions: getProperties, getStates, getProperty, setProperty and copy.

initialize() and run() functions

    void ModelExample::initialize(UByte d,State *s) {
        ConstitutiveModel::initialize(d,s);
        updateParameters();
    }

    void ModelExample::run(UByte d,State *s) {
        ConstitutiveModel::run(d,s);

        if (s->modulus_reduction_factor_ > 0.0)
            moduliReduction(s->modulus_reduction_factor_);
// excerpt-state-start
        if (s->state_ & shear_now) s->state_ |= shear_past;
        s->state_ &= ~shear_now;
        if (s->state_ & tension_now) s->state_ |= tension_past;
        s->state_ &= ~tension_now;
// excerpt-state-end
        UInt iPlas = 0;

        ModelElastic::elasticTrial(s);
        s->viscous_ = true;        
        if (!canFail()) return;

        SymTensorInfo info;
        DVect3 prin = s->stnS_.getEigenInfo(&info);

        Double fs  = nph_*prin.z() - prin.x() - csn_;
        Double fsd = fs / rc_;
        Double ftz = prin.z() - tension_;
        if (fsd > 0.0 && fsd >= ftz)
            shearCorrection(s,&prin,&iPlas,fs);
        else if (ftz > 0.0 && ftz >= fsd)
            tensionCorrection(s,&prin,&iPlas,ftz,brittle_);
        apexCorrection(friction_,s,&prin,&iPlas,brittle_);

        if (iPlas) {
            s->stnS_ = info.resolve(prin);
            s->viscous_ = false;
        }
    }

    bool ModelExample::updateParameters() {
        Double rsin = std::sin(friction_ * degrad);
        nph_ = (1.0 + rsin) / (1.0 - rsin);
        csn_ = 2.0 * cohesion_ * sqrt(nph_);
        rsin = std::sin(dilation_ * degrad);
        nps_ = (1.0 + rsin) / (1.0 - rsin);
        rc_  = std::sqrt(1.0 + nph_*nph_);
        //
        if (friction_) {
            Double apex = cohesion_ / std::tan(friction_ * degrad);
            tension_ = std::min(tension_,apex);
        }
        //
        ModelElastic::updateParameters();
        //
        Double ra = e1_ - nps_ * e2_;
        Double rb = e2_ - nps_ * e1_;
        Double rd = ra - rb * nph_;
        sc1_  = ra / rd;
        sc3_  = rb / rd;
        sc2_  = e2_ * (1.0 - nps_) / rd;
        //
        return true;
    }

    bool ModelExample::updateParameters(bool bEUpdated, Double *sf1, Double *sf3) {
        if (cohesion_ < 0.0)
            throw std::runtime_error("Example model: cohesion is not allowed less than 0.");

        if (friction_ >= 90.0 || friction_ < 0.0)
            throw std::runtime_error("Example model: friction angle is not in the valid range of 0 to 90.");

        if (dilation_ > friction_)
            throw std::runtime_error("Example model: dilationn angle is not allowed greater than friction angle.");

        Double rsin = std::sin(friction_ * degrad);
        nph_ = (1.0 + rsin) / (1.0 - rsin);
        csn_ = 2.0 * cohesion_ * sqrt(nph_);
        rsin = std::sin(dilation_ * degrad);
        nps_ = (1.0 + rsin) / (1.0 - rsin);
        rc_  = std::sqrt(1.0 + nph_*nph_);
        //
        if (friction_) {
            Double apex = cohesion_ / std::tan(friction_*degrad);
            tension_ = std::min(tension_,apex);
        }
        //
        if (!bEUpdated) ModelElastic::updateParameters();
        //
        Double ra = e1_ - nps_*e2_;
        Double rb = e2_ - nps_*e1_;
        Double rd = ra  - rb*nph_;
        sc1_ = ra/rd;
        sc3_ = rb/rd;
        sc2_ = e2_*(1.0 - nps_)/rd;
        //
        if (sf1) *sf1 = -1.0/rd;
        if (sf3) *sf3 = nps_/rd;
        //
        return !bEUpdated;
    }

FISH Support for Constitutive Models

The FISH intrinsic zone.prop is available in FLAC3D. This can be used on the left- or right-hand side of an expression. Thus,

val = zone.prop(zp, p_name)

stores in val the floating-point value of the property named p_name in the zone with address zp. p_name may be a string containing the property name, or a FISH variable that evaluates to such a string. For example, zone.prop(zp,'bulk') would refer to the bulk modulus. If there is no constitutive model in zp, or the model does not possess the named property, then 0.0 is returned. Similarly,

zone.prop(zp, p_name) = val

stores val in the property named p_name in zone zp. Nothing is stored if there is no constitutive model in zp, or if the model does not possess the named property, or val is not an integer or floating-point number. In both uses, zp must be a zone pointer and p_name must either be a string or a FISH variable containing a string.

Creating User-Written Model DLLs

Using the Template

The simplest way to create a constitutive model project in Visual Studio 2017 or 2019 is to use the ItascaProjectTemplates extension provided by Itasca.

Loading and Running User-Written Model DLLs

Model DLL files may be loaded into FLAC3D while it is running by giving the command zone cmodel load filename with the file name of the DLL. DLL files will be automatically loaded if they are placed in the “exe64\plugins\cmodel” folder. Thereafter, the new model name and property names will be recognized by FLAC3D and FISH functions that refer to the model and its properties. If the zone cmodel load command is given for a model that is already loaded, nothing will be done, but an informative message will be displayed.

Before constitutive model plugins can be assigned to zones, the model must be configured for their use by giving the model configure plugin command. Once configured, the model will not cycle unless your FLAC3D license includes the C++ plugin option.