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 “jmodels007.dll,” which defines the specific interface used by the joint 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 the following.

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 consistently-sized 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 × 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.

Example Constitutive Model

The source codes of all constitutive models used by 3DEC are provided for the user to inspect or adapt. Here are extracted, for illustration, parts of the Mohr-Coulomb elastic/plastic model contained in files “jmodelexample.*”. jmodelexample.h provides the class specification for the model, which also includes the definition of the model’s unique type number.

In this model, some of the variables are for internal use only (e.g., tan_friction_): they occupy memory in each subcontact, but they are not available for the user of 3DEC to change or print out. Also note that the getProperty()/setProperty interface is used for Save/Restore.

jmodelexample.h:

#pragma once

#include "jmodels/src/jointmodel.h"

namespace jmodels
{
  class JModelExample : public JointModel {
  public:
    JModelExample();
    virtual String         getName() const;
    virtual String         getPluginName() const { return getName(); }
    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 JModelExample *clone() const { return new JModelExample(); }
    virtual Double         getMaxNormalStiffness() const { return kn_; }
    virtual Double         getMaxShearStiffness() const { return ks_; }
    virtual void           copy(const JointModel *mod);
    virtual void           run(UByte dim,State *s); // If !isValid(dim) calls initialize(dim,s)
    virtual void           initialize(UByte dim,State *s); // calls setValid(dim)
    // Optional 
    virtual Double         getStressStrengthRatio(const Double &,const DVect3 &) const { return 10.0; }
    virtual void           scaleProperties(const Double &,const std::vector<UInt> &) { throw std::runtime_error("Does not support property scaling"); }
    virtual bool           supportsStressStrengthRatio() const { return false; }
    virtual bool           supportsPropertyScaling() const { return false; }
  private:
    Double kn_;
    Double ks_;
    Double cohesion_;
    Double friction_;
    Double dilation_;
    Double tension_;
    Double zero_dilation_;
    Double res_cohesion_;
    Double res_friction_;
    Double res_tension_;
    Double tan_friction_;
    Double tan_dilation_;
    Double tan_res_friction_;
    Int    kn_tab_;
    Int    ks_tab_;
  };
} // 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 (e.g., tan_friction_) are not computed at each cycle, rather during initialization in the Initialize method. Also note the use of the State structure in providing displacement increments and forces. Please refer to the file “jmodelexample.cpp” for listings of member functions: getProperties, getStates, getProperty, setProperty and copy.

`initialize()` and `run()` functions

  void JModelExample::initialize(UByte dim,State *s) 
  {
    JointModel::initialize(dim,s);
    tan_friction_    = tan(friction_ * dDegRad);
    tan_res_friction_ = tan(res_friction_ * dDegRad);
    tan_dilation_    = tan(dilation_ * dDegRad);
  }

  void JModelExample::run(UByte dim,State *s) 
  {
    JointModel::run(dim,s);
    /* --- state indicator:                                  */
    /*     store 'now' info. as 'past' and turn 'now' info off ---*/
    if (s->state_ & slip_now) s->state_ |= slip_past;
    s->state_ &= ~slip_now;
    if (s->state_ & tension_now) s->state_ |= tension_past;
    s->state_ &= ~tension_now;

    Double kna  = kn_ * s->area_;
    Double ksa  = ks_ * s->area_;

    // normal force
    Double fn0 = s->normal_force_;
    s->normal_force_inc_ = -kna * s->normal_disp_inc_;
    s->normal_force_ += s->normal_force_inc_;

    // correction for time step in which joint opens (or goes into tension)
    // s->dnop_ is part of s->normal_disp_inc_ at which separation or tension takes place
    s->dnop_ = s->normal_disp_inc_;
    if ((fn0 > 0.0)               && 
        (s->normal_force_ <= 0.0) && 
        (s->normal_force_inc_ < 0.0))
    {
        s->dnop_ = -s->normal_disp_inc_ * fn0 / s->normal_force_inc_;
        if (s->dnop_ > s->normal_disp_inc_) s->dnop_ = s->normal_disp_inc_;
    }

    // tensile strength
    Double ten;
    if (s->state_)
      ten = -res_tension_ * s->area_;
    else
      ten = -tension_ * s->area_;

    // check tensile failure
    bool tenflag = false;
    if (s->normal_force_ <= ten) 
    {
      s->normal_force_  = ten;
      if (!s->normal_force_)
      {
        s->shear_force_ = DVect3(0,0,0);
        tenflag = true; // complete tensile failure
      }
      s->state_ |= tension_now;
      s->normal_force_inc_ = 0.0;
      s->shear_force_inc_ = DVect3(0,0,0);
    }

    // shear force
    if (!tenflag) 
    {
      s->shear_force_inc_ = s->shear_disp_inc_ * -ksa;
      s->shear_force_ += s->shear_force_inc_;
      Double fsm = s->shear_force_.mag();
      // shear strength
      Double fsmax;
      if (!s->state_)
        fsmax = cohesion_ * s->area_ + tan_friction_ * s->normal_force_;
      else 
      { // if residual friction is zero, take peak value
        Double resamueff = tan_res_friction_;
        if (!resamueff) resamueff = tan_friction_;
        fsmax = res_cohesion_ * s->area_ + resamueff * s->normal_force_;
      }
      if (fsmax < 0.0) fsmax = 0.0;
      //  check for slip
      if (fsm >= fsmax) 
      {
        Double rat = 0.0;
        if (fsm) rat = fsmax / fsm;
        s->shear_force_ *= rat;
        s->state_ |= slip_now;
        s->shear_force_inc_ = DVect3(0,0,0);
        // dilation
        if (dilation_)
        {
          Double zdd = zero_dilation_;
          Double usm = s->shear_disp_.mag();
          if (!zdd) zdd = 1e20;
          if (usm < zdd) 
          {
            Double dusm  = s->shear_disp_inc_.mag();
            Double dil = 0.0;
            if (!s->state_) dil = tan_dilation_;
            else
            {
              // if residual dilation is zero, take peak value
              //     Double resdileff = tan_res_dilation_;
              // Note: In CLJ1 in 3DEC, no residual dilation is defined
              Double resdileff = tan_dilation_;
              if (!resdileff) resdileff = tan_dilation_;
              dil = resdileff;
            }
            s->normal_force_ += kna * dil * dusm;
          }
        } // dilation
      } // fsm>fsmax
    } // if (!tenflg)
  }

FISH Support for Constitutive Models

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

val = block.subcontact.prop(cxp, p_name)

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

block.subcontact.prop(cxp, p_name) = val

stores val in the property named p_name in subcontact cxp. Nothing is stored if there is no constitutive model in cxp, or if the model does not possess the named property, or val is not an integer or floating-point number. In both uses, cxp must be a subcontact 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.

To see how it can be done manually, see Manual Configurations in the |flac3d| Implementation document.

Loading and Running User-Written Model DLLs

Joint model DLL files may be loaded into 3DEC while it is running by giving the command contact jmodel load filename with the file name of the DLL. DLL files will be automatically loaded if they are placed in the “exe64\plugins\jmodel” folder. Thereafter, the new model name and property names will be recognized by 3DEC and FISH functions that refer to the model and its properties. If the contact jmodel 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 the current 3DEC license includes the C++ plugin option.

⇄

⇐ Methodology | Programmer’s Interface ⇒

Was this helpful? ...

Updated: Feb 25, 2024