Implementation

Utility Structures

A few structures/classes are provided to help in writing and communicating with constitutive models. Some are provided by “base009.dll”, which defines the base level of functionality common to all plugin interfaces. Others are provided by “conmodel009.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.

“base009.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 assigned 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 “conmodel009.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 "conmodel.h"
#include "convert.h"

namespace models {
    class ModelExample : public ConstitutiveModel {
      public:
        ModelExample(unsigned short option = 0);
        string getName() const override;
        string getFullName() const override;
        uint32 getMinorVersion() const override;
        string getProperties() const override;
        string getStates() const override;
        base::Property getProperty(uint32 index) const override;
        void setProperty(uint32 index, const base::Property& p, uint32 restoreVersion = 0) override;
        ConstitutiveModel* clone() const override { return new ModelExample(); }
        double getConfinedModulus() const override { return bulk_ + shear_*d4d3; }
        double getShearModulus() const override { return shear_; }
        double getBulkModulus() const override { return bulk_; }
        void copy(const ConstitutiveModel* mod) override;
        void run(uint32 dim, State* s) override;
        void initialize(uint32 dim, State* s) override;
        // Optional
        bool isPropertyAdvanced(uint32 i) const override;
        bool supportsStrengthStressRatio() const override { return true; }
        bool supportsPropertyScaling() const override { return true; }
        double getStrengthStressRatio(const SymTensor& st) const override;
        void scaleProperties(const double& scale, const std::vector<uint32>& props) override;
        bool isKGEv() const override { return false; }

      private:
        double moduliReduction(const double& factor);
        void apexCorrection(const double& fric, State* s, DVect3* prin, uint32* iPlasticity = nullptr, bool bBrittle = false);
        void tensionCorrection(State* s, DVect3* prin, uint32* iPlasticity, const double& ftz, bool bBrittle = false);
        void shearCorrection(State* s, DVect3* prin, uint32* iPlasticity, const double& fs);

        double bulk_ = 0.0;
        double shear_ = 0.0;
        double cohesion_ = 0.0;
        double friction_ = 0.0;
        double dilation_ = 0.0;
        double tension_ = 0.0;
        bool brittle_ = false;

        double e1_ = 0.0;
        double e2_ = 0.0;
        double g2_ = 0.0;
        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 initialization 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(uint32 d,State *s) {
        ConstitutiveModel::initialize(d,s);

        e1_ = bulk_ + shear_ * d4d3;
        e2_ = bulk_ - shear_ * d2d3;
        g2_ = 2.0 * shear_;

        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);
        }
        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;
    }

    void ModelExample::run(uint32 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		
        uint32 iPlas = 0;

        double e11 = s->stnE_.s11();
        double e22 = s->stnE_.s22();
        double e33 = s->stnE_.s33();
        s->stnS_.rs11() += e11 * e1_ + (e22 + e33) * e2_;
        s->stnS_.rs22() += (e11 + e33) * e2_ + e22 * e1_;
        s->stnS_.rs33() += (e11 + e22) * e2_ + e33 * e1_;
        s->stnS_.rs12() += s->stnE_.s12() * g2_;
        s->stnS_.rs13() += s->stnE_.s13() * g2_;
        s->stnS_.rs23() += s->stnE_.s23() * g2_;

        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;
        }
    }

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.

Loading and Running User-Written Model DLLs

After the model DLL file is created (Writing New Constitutive Models), it 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 plugins command. Once configured, the model will not cycle unless your FLAC3D license includes the C++ plugin option.