Implementation

Utility Structures

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

“base005.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<char>, 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 “conmodel006.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 private variables than property names (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 "../src/conmodel.h"

namespace models {
    class ModelExample : public ConstitutiveModel {
    public:
        ModelExample();

        virtual String        getName() const;
        virtual String        getFullName() const;
        virtual UInt          getMinorVersion() const; 
        virtual String        getProperties() const;
        virtual String        getStates() const;
        virtual Variant       getProperty(UInt index) const;
        virtual void          setProperty(UInt index,const Variant &p,UInt restoreVersion=0);
        virtual ModelExample *clone() const { return new ModelExample(); }
        virtual Double        getConfinedModulus() const { return bulk_+shear_*d4d3; }
        virtual Double        getShearModulus() const { return shear_; }
        virtual Double        getBulkModulus() const { return bulk_; }
        virtual bool          supportsHystereticDamping() const { return true; }
        virtual void          copy(const ConstitutiveModel *mod);
        virtual void          run(UByte dim,State *s); 
        virtual void          initialize(UByte dim,State *s);
        // Optional
        virtual bool          isPropertyAdvanced(UInt i) const;
        virtual Double        getStressStrengthRatio(const SymTensor &st) const;
        virtual void          scaleProperties(const Double &scale,const std::vector<UInt> &props);
        virtual bool          supportsStressStrengthRatio() const { return true; }
        virtual bool          supportsPropertyScaling() const { return true; }
    private:
        Double bulk_,shear_,cohesion_,friction_,dilation_,tension_;
        Double e1_,e2_,g2_,nph_,csn_,sc1_,sc2_,sc3_,bisc_,e21_,rnps_;
        bool brittle_;
    };
} // namespace models

// EOF

Constant Definitions

// Plasticity Indicators
static const UInt shear_now    = 0x01;
static const UInt tension_now  = 0x02;
static const UInt shear_past   = 0x04;
static const UInt tension_past = 0x08;

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, private model variables e1_, e2_, g2_, etc. are not computed at each cycle. Also note the use of the State structure in providing strain increments and stresses. In general, separate sections should 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 dim,State *s)  {
    ConstitutiveModel::initialize(dim,s);
    e1_ = bulk_ + shear_*d4d3;
    e2_ = bulk_ - shear_*d2d3;
    g2_ = shear_*2.0;
    Double rsin = std::sin(friction_ * degrad);
    nph_  = (1.0 + rsin) / (1.0 - rsin);
    csn_  = 2.0 * cohesion_ * sqrt(nph_);
    if (friction_) {
        Double apex = cohesion_ / std::tan(friction_ * degrad);
        tension_ = std::min(tension_,apex);
    }
    rsin = std::sin(dilation_ * degrad);
    rnps_ = (1.0 + rsin) / (1.0 - rsin);
    Double ra = e1_ - rnps_ * e2_;
    Double rb = e2_ - rnps_ * e1_;
    Double rd = ra - rb * nph_;
    sc1_  = ra / rd;
    sc3_  = rb / rd;
    sc2_  = e2_ * (1.0 - rnps_) / rd;
    bisc_ = std::sqrt(1.0 + nph_*nph_) + nph_;
    e21_  = e2_ / e1_;
}

void ModelExample::run(UByte dim,State *s) {
    ConstitutiveModel::run(dim,s);
    if (s->modulus_reduction_factor_>0.0) {
        Double shear_new = shear_ * s->modulus_reduction_factor_;
        e1_ = bulk_ + shear_new * d4d3;
        e2_ = bulk_ - shear_new * d2d3;
        g2_ = 2.0 * shear_new;
        Double ra = e1_ - rnps_ * e2_;
        Double rb = e2_ - rnps_ * e1_;
        Double rd = ra - rb * nph_;
        sc1_  = ra / rd;
        sc3_  = rb / rd;
        sc2_  = e2_ * (1.0 - rnps_) / rd;
        e21_  = e2_ / e1_;
    }
    // plasticity indicator:
    // store 'now' info. as 'past' and turn 'now' info off
    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;
    UInt plas = 0;

    /* --- trial elastic stresses --- */
    Double e11 = s->stnE_.s11();//strain tensors normal components
    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_;

    // default settings, altered below if found to be failing
    s->viscous_ = true;  // Allow stiffness-damping terms

    if (canFail()) {
        // Calculate principal stresses
        SymTensorInfo info;
        DVect3 prin = s->stnS_.getEigenInfo(&info);

        /* --- Mohr-Coulomb failure criterion --- */
        Double fsurf = prin.x() - nph_ * prin.z() + csn_;//
        /* --- Tensile failure criteria --- */
        Double tsurf = tension_ - prin.z();//
        Double pdiv = -tsurf + (prin.x() - nph_ * tension_ + csn_) * bisc_;

        /* --- tests for failure */
        if (fsurf < 0.0 && pdiv < 0.0) {
            plas = 1;
            /* --- shear failure: correction to principal stresses ---*/
            s->state_ |= shear_now;
            prin.rx() -= fsurf * sc1_;
            prin.ry() -= fsurf * sc2_;
            prin.rz() -= fsurf * sc3_;
        } else if (tsurf < 0.0 && pdiv > 0.0) {
            plas = 2;
            /* --- tension failure: correction to principal stresses ---*/
            s->state_ |= tension_now;
            Double tco = e21_ * tsurf;
            prin.rx() += tco;
            prin.ry() += tco;
            prin.rz()  = tension_;
        }

        if (plas) {
            s->stnS_ = info.resolve(prin);// transform back to refrence frame
            s->viscous_ = false; // Inhibit stiffness-damping terms
        }
    }
}

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 2010 is to use the Project Template add-in provided by Itasca.

Manual Configurations

In order to create a DLL in Visual Studio 2010, it is first necessary to create a solution. The solution will contain projects that are essentially a collection of C++ source and header files and their dependencies.

An example solution and project have been provided. This collection of files can be found in the “pluginfiles\models\example” folder of the FLAC3D application data directory (“My Documents\Itasca\FLAC3D600” by default). To create a custom DLL, the user may simply wish to modify the example provided. We do not, however, recommend this, since the original example will no longer be available for reference. Additionally, Visual Studio embeds a unique identifier (a GUID) in each solution and project file. Copying a project and renaming it can cause serious confusion within the IDE. Instead, we recommend you create a new project using the following steps and settings.

1. Launch Visual Studio 2010 from the Start menu and select File ‣ New ‣ Project. Under “Project Types,” select “Visual C++/Win32”. Under “Templates,” select “Win32 project”. Select a location and name for the project, and press OK.
2. A Win32 Application Wizard dialog will appear. Select “Application Setting” on the left. Under “Application Type,” select “DLL”. Under “Additional Options,” select “Empty Project”. Then press Finish.
3. Copy the files “modelexample.h”, “modelexample.cpp”, “version.rc” and “version.txt” from “pluginfiles\models\examples” into the directory created by Visual Studio for your solution and project. Rename the C++ source files to indicate your constitutive model (for example, “modeltest.h” and “modeltest.cpp”).
4. In the Solution Explorer window of Visual Studio, right-click on the “Header Files” folder and add “modeltest.h”. Then right-click on the “Source Files” folder and add “modeltest.cpp”. Add the file “version.rc” to the “Resources” folder. Add the file “version.txt” as a project object. (“Add Existing Item” from the project context menu.) These files can now be edited to create your custom constitutive model.
5. Right-click on the project entry and select “Properties”. Make certain that “Configuration:” reads “Active(Debug)” and “Platform:” reads “Active(x64)”. You can create Release and x64 versions of the project later by copying these steps with appropriate modifications.
6. Under “Configuration Properties,” select the “General” entry. Change the “Output Directory” entry to “$(ProjectDir)$(ConfigurationDebug)”.
7. Under “C/C++,” select the “General” entry. Add the directories “interface” and “cmodels\src” (separated by ;), both sub-directories of the FLAC3D install directory, to “Additional Include Directories”. In a typical install, this entry might read “C:\Program Files\Itasca\FLAC3D600\pluginfiles\interface”; “C:\Program Files\Itasca\FLAC3D600\pluginfiles\cmodels\src”. You should also change “Warning Level” to Level 4, and set “Treat Warnings As Errors” to Yes.
8. Under “C/C++,” select the “Preprocessor” entry. Change “Preprocessor Definitions” so that “_DEBUG” now reads “MODELDEBUG.”
9. Under “C/C++,” select the “Code Generation” entry. Change “Runtime Library” to “Multi-threaded DLL (/MD).”
10. Select “Linker/General”. The name of the DLL produced should follow the convention “model<name>006_64.dll”, where <name> is the (unique) name of your constitutive model as returned by the getName() method. The “model” prefix identifies it as a constitutive model, and 006 indicates the major version number, and therefore binary compatibility. For example, this entry should read “$(OutDir)\modeltestd006_64.dll”. The “d” is to allow you to distinguish a debug version of the model from a release version.
11. Select “Linker/Input”. The “Additional Dependencies” should add the files “base005_64.lib” and “conmodel006_64.lib”, separated by a space, and including the FLAC3D install directory. For instance, in a typical install this would read “C:\Program Files\Itasca\FLAC3D600\pluginfiles\lib\exe64\base005_64.lib” and “C:\Program Files\Itasca\FLAC3D600\pluginfiles\lib\exe64\conmodel006_64.lib”.

To create a release version, make the same property changes in the “Release” configuration properties, but do not add the “MODELDEBUG” preprocessor definition. Also, remove the “d” suffix from the model name.

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.