Difference between revisions of "Creating a Custom Imperfect Interface Law"

From OSUPDOCS
Jump to navigation Jump to search
 
(17 intermediate revisions by the same user not shown)
Line 1: Line 1:
This section explains how to write <tt>C++</tt> code create a custom [[Imperfect Interfaces|imperfect interface]] for use in [[NairnMPM]] along with available [[Contact Laws#Imperfect Interface Contact Laws|imperfect interface contact laws]]..
This section explains how to write <tt>C++</tt> code create a custom [[Imperfect Interfaces|imperfect interface]] for use in [[NairnMPM]] along with available [[Contact Laws#Imperfect Interface Contact Laws|imperfect interface contact laws]].


== Getting Started ==
== Getting Started ==
Line 7: Line 7:
=== Class Source Code ===
=== Class Source Code ===


Custom imperfect interface should start by duplicating <tt>NonlinearInterface.cpp</tt> and <tt>NonlinearInterface.hpp</tt> files (or the files for one of its subclasses). All custom lows should inherit from the <tt>NonlinearInterface</tt> class. The new contact law will be used along with other [[Contact Laws]] and therefore will need a unique name and ID (an integer). In the new contact law's header file, replace <tt>NONLINEARINTERFACELAW</tt> copies above with a defined constant representing the new contact law's ID (which by convention is in UPPERCASE) and replace the number in the new constant's definition with the new ID. For example:
Custom imperfect interface should start by duplicating <tt>NonlinearInterface.cpp</tt> and <tt>NonlinearInterface.hpp</tt> files (or the files for one of its subclasses). All custom lows should inherit from the [[Nonlinear Imperfect Interface|<tt>NonlinearInterface</tt> class]]. The new contact law will be used along with other [[Contact Laws]] and therefore will need a unique name and ID (an integer). In the new contact law's header file, replace <tt>NONLINEARINTERFACELAW</tt> copies above with a defined constant representing the new contact law's ID (which by convention is in UPPERCASE) and replace the number in the new constant's definition with the new ID. For example:


  #define MYINTERFACE 201
  #define MYINTERFACE 201
Line 50: Line 50:
</ol>
</ol>


== Writing the Material Class Source Code ==
== Writing the Contact Law Class Source Code ==


The remaining steps can usually all be done by editing only the material's source code file. The requirements and optional methods are described in the following sections. This stage involves editing the template file created above and customizing the listed methods and/or deleting optional methods that are not needed. All methods that are overridden should be declared as <tt>virtual</tt> in the header file to insure the correct code is used at run time.
The remaining steps can usually all be done by editing only the contact law's source code file. The requirements and optional methods are described in the following sections. This stage usually starts by keeping by deleting all methods in the copied file that are not listed below and then changing the class name in each definition (''e.g''., change <tt>NonlinearInterface::</tt> to <tt>MyInterface::</tt> in each method). Any methods deleted from the <tt>cpp</tt> file should have its template deleted from the <tt>hpp</tt> file as well.


=== Creating Material Properties ===
=== Creating Material Properties ===


Usually a newly-created material type will have additional material properties. To create such properties, you need to define them, allow them in input command files (''i.e.'', update the <tt>DTD</tt> file), set them in the material class file, validate them, and finally use them in calculations. The following steps are needed:
Usually a newly-created contact law will have additional material properties than the parent law's (although reuse parent law properties when possible). To create such properties, you need to define them, allow them in input command files (''i.e.'', update the <tt>DTD</tt> file), set them in the contact law's class file, and use them in calculations. The following steps are needed:


<ol>
<ol>
Line 72: Line 72:
</ol>
</ol>


Once a property variable and name are defined, you set that variable, check it setting, and use it with the following methods:
One properties are signed, they are initialized and printed with the following methods:


==== Required Initialization Methods ====
==== Required Initialization Methods ====


; <tt>char *InputMaterialProperty(char *xName,int &input,double &gScaling)</tt>
; <tt>MyInterface::MyInterface(char *matName,int matID) : NonlinearInterface(matName,matID)</tt>
: If <tt>xName</tt> string matches a property name for this material, set <tt>input</tt> to <tt>DOUBLE_NUM</tt> or <tt>INT_NUM</tt> (depending on the type of variable) and return a pointer to the class variable for that property. If <tt>gScaling</tt> is set, the input  value will be multiplied by <tt>gScaling</tt> after it is read. If <tt>xName</tt> is not a recognized property return <tt>Parent::InputMat(xName,input)</tt> to allow the super class to look for their valid property types (replace <tt>Parent</tt> with the name of the superclass).
: This method is the default constructor. It can initialize any custom interface properties to default values. It should als call constructor for the parent contact law class. If no defaults are needed, this method can be omitted. The <tt>MyInterface:::</tt> is the custom class name. This term is omitted from subsequent methods names (but is needed).


; <tt>const char *VerifyAndLoadProperties(int np)</tt> (optional)
; <tt>char *InputContactProperty(char *xName,int &input,double &gScaling)</tt>
: This method is called after the input file is read but before the new material is printed to the results file. You can use this method to verify the input material properties are physically allowed. If not, return a string with an error message and the simulation will be aborted. If there are no errors, return <tt>Parent::VerifyAndLoadProperties(np)</tt> to let super class check properties as well (replace <tt>Parent</tt> with the name of the superclass) The input <tt>np</tt> is a constant for the analysis type (''e.g.'', plane stress, plane strain, 3D, ''etc.''). If the properties are valid, but maybe not allowed in current MPM mode, that check should be done instead in <tt>ValidateForUse()</tt> below. This method is only called once at the beginning. For efficient, it is therefore a good place to to calculate material properties that are independent of the particle state and thus will remain constant throughout the calculation (''e.g.'', to calculate specific properties by dividing by density or to convert to units more convenient to the constitutive law).
: If <tt>xName</tt> string matches a property name for this material, set <tt>input</tt> to <tt>DOUBLE_NUM</tt> or <tt>INT_NUM</tt> (depending on the type of variable) and return a pointer to the class variable for that property. If <tt>gScaling</tt> is set, the input  value will be multiplied by <tt>gScaling</tt> after it is read. If <tt>xName</tt> is not a recognized property return <tt>NonlinearInterface::InputMat(xName,input)</tt> to allow the super class to look for their valid property types (replace <tt NonlinearInterface</tt> with different name if you subclassed some other contact law).


; <tt>void PrintMechanicalProperties(void)</tt>
; <tt>void PrintContactLaw(void) const</tt>
: This method should print all mechanical properties or call a super class and thenprint just the new mechanical properties. Use a style similar to other materials. To help in formatting, you can use the <tt>MaterialBase</tt> class method <tt>PrintProperty(label,prop,units)</tt> to print a label, a property numeric value, and units within one column. You can use several calls in sequence to get several properties on the same line. You can also call <tt>PrintProperty(text,right)</tt> to print text in a column and right or left justified if <tt>right</tt> is true or false. This method need not pass on to <tt>MaterialBase</tt> because it does not print any mechanical properties. This method is called after <tt>VerifyAndLoadProperties()</tt> is done.
: This method should print all properties or call a super class and then print just the new law properties. Use a style similar to other contact laws. To help in formatting, you can use the <tt>MaterialBase</tt> class method <tt>PrintProperty(label,prop,units)</tt> to print a label, a property numeric value, and units within one column. You can use several calls in sequence to get several properties on the same line. You can also call <tt>PrintProperty(text,right)</tt> to print text in a column and right or left justified if <tt>right</tt> is true or false.
 
==== Optional Initialization Methods ====
 
; <tt>void ValidateForUse(int np)</tt>
: This method is called just before the first MPM time step and only for materials used by one or more materials. Throw a <tt>CommonException</tt> if this material cannot by used in current MPM mode (type specified in <tt>np</tt> which will be <tt>PLANE_STRAIN_MPM</tt>, <tt>PLANE_STRESS_MPM</tt>, <tt>AXISYMMETRIC_MPM</tt>, or <tt>THREED_MPM</tt>). If no exceptions, must call the same method in the parent class. Basic material properties are usually checked in <tt>VerifyAndLoadProperties()</tt> instead; this method is for materials that may have valid properties, but may be contingent on other MPM settings or for materials only implement for certain styles of simulations.
 
; <tt>FillTransportProperties(TransportProperties *)</tt>
: This method is called by <tt>MaterialBase::VerifyAndLoadProperties()</tt> and can be used calculate transport properties that are independent of the particle state and thus will remain constant throughout the calculation. The <tt>MaterialBase</tt> class automatically finds diffusion and conductivity tensors for isotropic materials by using material properties <tt>D</tt> (in <tt>diffusionCon</tt>) for diffusion constant and <tt>kCond</tt> (in <tt>kCond</tt>) for thermal conductivity. It only needs to be overridden if something more is needed.
 
; <tt>void PrintTransportProperties(void)</tt>
: This method should print all transport properties in format similar to other materials (see help <tt>PrintProperty()</tt> methods described above). It should only print them if transport is activated (''i.e.'', <tt>if(DiffusionTask::active)</tt> or <tt>if(ConductionTask::active)</tt>). The <tt>MaterialBase</tt> class prints isotropic properties and the <tt>ORTHO</tt> and <tt>TRANSISO1(2)</tt> classes print anisotropic properties. No additional printing is needed if one of these classes handles the task.
 
; <tt>void SetInitialParticleState(MPMBase *mptr,int np,int offset) const</tt>
; This method is called in preliminary calculations before a simulation states. It can be used, if needed, to set and material properties for each material point (in the <tt>mptr</tt> pointer).
 
=== Basic Class Accessors ===
 
These methods (required unless specified as optional) return basic facts about the material:


==== Required Accessors ====
==== Required Accessors ====


; <tt>const char *MaterialType(void)</tt>
; <tt>const char *MaterialType(void)</tt>
: Return a short name to describe the material. This string is printed with material properties in the output of simulations.
: Return a short name to describe the imperfect interface law. This string is printed with contact law properties in the output of simulations.
 
; <tt>double WaveSpeed(bool,MPMBase *)</tt>
: Return the maximum wave speed for material in [[ConsistentUnits Command#Legacy and Consistent Units|velocity units]]. This method is called once for each material point at start of calculation and after material properties have been defined. If the wave speed might change during the simulation, be conservative and return the maximum possible save speed. The first parameter is true for 3D calculations or false for 2D. This speed is needed to pick the time step for explicit calculations. The second parameter is pointer to a material point in case the wave speed depends on the initial particle state. This method is also used in some crack propagation methods. If the material supports crack propagation, make sure it gives good results for these uses.
 
==== Optional Accessors ====
 
; <tt>double CurrentWaveSpeed(bool,MPMBase *)</tt>
: Return the current wave speed, which might depend on particle state. This method is only used by the [[AdjustTimeStep Custom Task|<tt>AdjustTimeStep</tt> task]] to change time step as needed; without this method, the time step will not adjust for this new material.
 
; <tt>double ShearWaveSpeed(bool)</tt>
: Return the shear wave speed for the material. This method is only used by [[Boundary Condition Styles#Silent Boundary Conditions|silent boundary conditions]] for force. The <tt>MaterialBase</tt> class returns <tt>WaveSpeed()/sqrt(3)</tt>. You only need to override this method if you want a better value. Also note that [[Boundary Condition Styles#Silent Boundary Conditions|silent boundary conditions]]  only work for isotropic materials.
 
; <tt>double MaximumDiffusion(void)</tt>
: Return maximum diffusion constant in [[ConsistentUnits Command#Legacy and Consistent Units|diffusion units]]. This method is called once for each material at the start of the calculations and after material properties are defined. The <tt>MaterialBase</tt> class returns the appropriate result for isotropic materials. You only need to override if you need a different result.
 
; <tt>double MaximumDiffusivity(void)</tt>
: Return maximum thermal diffusivity in [[ConsistentUnits Command#Legacy and Consistent Units|diffusion units]] = k/(100 rho Cp). This method is called once for each material at the start of the calculations and after material properties are defined. The <tt>MaterialBase</tt> class returns the appropriate result for isotropic materials. You only need to override if you need a different result.
 
; <tt>double GetCurrentRelativeVolume(MPMBase *)</tt>
: Return current relative volume to be used to convert tracked stress to true stress. Low strain materials track stress/initial density and they should return 1 (which is what base material class returns). Large deformation materials should be tracking specific Cauchy stress (= Kirchoff stress/initial density) and therefore should return relative volume (''i.e.'', initial density/current density).
 
; <tt>Tensor GetStress(Tensor *sp,double pressure)</tt>
: Return current stress. Most materials just return the contents of <tt>sp</tt> which is the particle stress. Materials that track pressure and deviatoric stress separately (or use some other stress tracking scheme), should reconstruct and return the full stress in this method.
 
; <tt>bool SupportsArtificialViscosity(void)</tt>
: Materials that implement [[Common Material Properties#Artificial Viscosity|artificial viscosity]] must override this method to return <tt>TRUE</tt> and then implement the calculations in the constitutive law methods.
 
; <tt>bool SupportsDiffusion(void) const</tt>
: if a material does not support diffusion or poroelasticity, override and return false only called when diffusion is active. The <tt>MaterialBase</tt> class returns true.
 
; <tt>double GetMaterialConcentrationSaturation(MPMBase *mptr) const</tt>
: Return the saturation concentration for the material. Only needed if the value depends on state of the material point. The <tt>MaterialBase</tt> returns the saturation concentration property.
 
; <tt>double GetRho(MPMBase *) const</tt>
: Return the density for the material. Only needed if the value depends on state of the material point. The <tt>MaterialBase</tt> returns the density property.
 
=== Alternate Particle Strain ===
 
The particle information can store two strains - total strain and an alternate strain. A given material type can use the alternate strain anyway it needs, but only certain styles are fully supported by visualization tools. The following table lists supported alternate strains and how visualization options plots elastic and plastic strains.
 
{| class="wikitable"
|-
! Type !! Stored Alternate Strain !! Viz Elastic Strain || Viz Plastic/Alternate Strain
|-
| <tt>NOTHING</tt> || nothing || total strain || 0
|-
| <tt>ENG_BIOT_PLASTIC_STRAIN</tt> || plastic strain || total - plastic strain || plastic strain
|-
| <tt>LEFT_CAUCHY_ELASTIC_B_STRAIN</tt> || elastic '''B''' tensor || elastic from '''B''' || total - elastic from '''B'''
|-
| <tt>LEFT_CAUCHY_TOTAL_B_STRAIN</tt> || total '''B''' tensor || total strain || 0
|-
| <tt>MEMBRANE_DEFORMATION</tt> || membrane deformation || total strain || 0
|}
 
Whenever strain is archived, total rotational strain is also archived. Thus all materials can plot total strain and any components of the total deformation gradient. Note that versions of [[NairnMPM]] prior to 11.0 used to have option for materials to partition elastic and plastic strains into the two particle strains. Such materials should only be visualized in the [[NairnFEAMPMViz]] tool provided with that old version. The use of current visualization tools on output from old calculations will not correctly plot elastic and plastic strains.
 
When creating a new material, the following method should return how it uses the alternate strain:
 
; <tt>int AltStrainContains(void)</tt> (optional)
: If the alternate strain tensor is used, return the quantity stored by the material. The supported options are in the "Type" column in the above table. If you choose to use that alternate strain for some other purpose, the material can return <tt>NOTHING</tt> to avoid confusing visualization calculations.
 
This above method is only needed in calculations to support correct [[MPM Global Archiving Options|global archiving]] and [[VTKArchive Custom Task|VTK archiving]] of strain quantities. The separate [[NairnFEAMPM]] and [[NairnFEAMPMViz]] tools evaluate strains from raw total and alternate strains depending on material type. A problem arises when adding your own materials. These tools will not know what is in the plastic strain. The current tools assume all unknown materials store <tt>ENG_BIOT_PLASTIC_STRAIN</tt> in the alternate strain. You can visualize results as follows:
 
# To see what is in particle strain, plot total strain.
# To see what is in alternate strain, plot plastic strain.
# Plots of elastic strain (which are difference of the two strains) may not make sense.
 
=== History Dependent Properties ===


For some materials, the constitutive law will depend on the state of the particle in the form of history-dependent data (e.g. cumulative plastic strain for plasticity or strain history for viscoelasticity). Such properties are implemented in a material by using history data. The methods in the material class are as follows:
=== MPM Step Calculations ===


; <tt>char *InitHistoryData(char *pchr,MPMBase *mptr)</tt>
The remaining source code is to use the law in MPM calculations. Much of the work is done in the parent <tt>NonlinearInterface</tt> class. Your custom subclass only needs to return details about the contact law it self. The current implementation assumes a decoupled law that the normal and tangential forces are given by:
: Tthis method is called once for each material point using this material at the start of the calculations. If <tt>pchr==NULL</tt>, It should allocate memory to store history data on the particle, otherwise assume <tt>pchr</tt> points to enough memory for history data. Fill history data and return pointer to the data. It can be as simple as a single variable or a pointer to an array of any number of variables.


; <tt>double GetHistory(int num,char *ptr)</tt>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
: This method should extract history variable number num from the allocated data pointer in <tt>ptr</tt>. The <tt>MaterialBase<tt> class assumes history is an array of doubles of length equal to <tt>NumberOfHistoryDoubles()</tt>.
<math>T_n = f_n([u_n]) \qquad {\rm and } \qquad T_t = f_t([u_t])</math>


; <tt>int NumberOfHistoryDoubles(void) const</tt>
or where normal traction <math>T_n</math> is a function only of normal interfacial discontinuity <math>[u_n]</math> and tangential traction <math>T_t</math> is a function only of tangential interfacial discontinuity <math>[u_t]</math>. This limitation may change in the future. As a result, each method defined below has two version - one for normal tractions and one for tangential tractions.
: Return number of doubles used in a simple array of doubles.


== MPM Step Calculations ==
==== Check for Stability ====


The bulk of MPM calculations for a material will be done in the material-dependent code that is called during each MPM time step. Below is a summary of the key methods. These methods should be made as efficient as possible.
; <tt>int CheckDnStability(double d,double m,double delta) const</tt><br><tt>int CheckDtStability(double d,double m,double delta) const</tt>
 
=== Prepare for Updates ===
 
These methods are called prior to constitutive law methods to allow update of material properties when those material properties depend on the current state of the particle (e.g., properties that depend on temperature, moisture, orientation, or current particle's stress or stain) and to support parallel code for strain updates.
 
Need to explain particle-dependent variables
 
; <tt>int SizeOfMechanicalProperties(int &altBufferSize) const</tt>
: The size
: The size


; <tt>void *GetCopyOfMechanicalProps(MPMBase *mptr,int np,void *matBuffer,void *altBuffer) const</tt>
==== Interface Law Properties ====
: This method is called just before the constitutive law on each time step. You can set any parameters for the material that depend on the current state of the particle. Things that never change (''i.e.'', properties that are independent of particle state) should be calculated in <tt>VerifyAndLoadProperties()<tt> instead. Some materials also used in FEA make use of <tt>LoadMechProps()</tt> in this task, but that method is limited to 2D because it only inputs a single rotation about the z axis.
 
; <tt>void GetTransportProps(MPMBase *,int,TransportProperties *)</tt> const (optional)
: This method is called when looping over material points to store parameters needed in transport calculations (''i.e.'', diffusion and conductivity tensors). It is called prior to the transport task to <tt>AddForces()</tt>. It is only needed for anistropic materials or for materials whose transport properties change depending on particle state. This method should load the required values into <tt>diffusionTensor</tt> and <tt>kCondTensor</tt> variables. It is automatically handled for subclasses of <tt>TRANSISO1(2)</tt> if the only effect is the current orientation. Transport properties that never change (''i.e.'', independent of particle state) should be calculated in <tt>FillTransportProperties()</tt> instead.
 
; <tt>double GetHeatCapacity(MPMBase *)</tt> (optional)
: Return current constant-strain heat capacity in the material (and code should always call this method when heat capacity is needed rather then directly using the heat capacity variable). This call allows materials to implement a state-dependent heat capacity. Return current heat capacity [[ConsistentUnits Command#Legacy and Consistent Units|heat capacity units]]  or in of nJ/(g-K) of using Legacy units. This method can be omitted if heat capacity is a constant.
 
; <tt>double GetCpMinusCv(MPMBase *mptr)</tt> (optional)
: Return the difference between constant-stress ( C<sub>p</sub>) and constant-strain ( C<sub>v</sub>) heat capacities (and code should always call this method when heat capacity is needed rather then directly using the heat capacity variable). Return the difference in[[ConsistentUnits Command#Legacy and Consistent Units|heat capacity units]]  or in of nJ/(g-K) of using Legacy units. This call is only used by conduction calculations because heat flow is based on C<sub>p</sub> instead of on C<sub>v</sub>. For an elastic material, this difference is given by -'''M'''.'''&alpha;'''T/&rho;, where '''M''' is stress-temperature tensor and '''&alpha;''' is thermal expansion tensor. For an isotropic material, this reduces to 3K&alpha;<sup>2</sup>T/&rho; where K is bulk modulus and &alpha; is linear thermal expansion coefficient. For an ideal gas, it reduces to nR/&rho;. If this method is omitted, the difference defaults to zero (which is often close enough for most solid materials).
 
=== Constitutive Law Methods ===


Once the properties are set (from previous section), the most important method is the one that implements the constitutive law:
; <tt>double GetFn(double delta,double area) const</tt><br><tt>double GetFt(double delta,double area) const</tt>
: These methods (one for normal and one for tangential directions) should return the force due to the traction law at the provided the displacement discontinuity <tt>delta</tt>. The force must be in [[ConsistentUnits Command#Legacy and Consistent Units|Force units]] for given displacement in [[ConsistentUnits Command#Legacy and Consistent Units|Length units]]. Many imperfect interface laws give a traction or stress. To convert such a law to force, multiply the traction by the contact area (provided in the <tt>area</tt>).


; <tt>void MPMConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np)</tt>
; <tt>double GetFnPrime(double delta,double area) const</tt><br><tt>double GetFtPrime(double delta,double area) const</tt>
: This method applies constitutive law for the material and updates all needed particle properties. The required updates include stress (should be a specific stress), strain, plastic strain, rotational strain, work energy, residual energy, plastic energy, heat energy and any history dependent variables defined for the material. To support thermal and solvent expansion, include their effect on the constitutive law. The input '''<tt>du</tt>''' is the velocity gradient times the time step, which gives displacement gradient for this time step. See note below on particle temperature updates.
: These methods (one for normal and one for tangential directions) should return derivative of the force due to the traction law at the provided the displacement discontinuity <tt>delta</tt>. The derivative must be in [[ConsistentUnits Command#Legacy and Consistent Units|Force units/Length units]] for given displacement in [[ConsistentUnits Command#Legacy and Consistent Units|Length units]]. Many imperfect interface laws give a traction or stress which have derivate for stress per unit length. To convert such a law proper derivative units, multiply the traction derivative by the contact area (provided in the <tt>area</tt>).


The update should never change particle temperature, because that would invalidate some thermodynamics calculations in [[NairnMPM]]. To cause a particle temperature change, call <tt>IncrementHeatEnergy()</tt> with a temperature or energy increment. The code will take care of translating this result into particle temperature rise when appropriate. The particle temperature will rise during adiabatic calculations. For isothermal calculations, the temperature will not rise, but the corresponding energy that was released will be reflected in heat energy. See [[Thermal Calculations#Thermodynamics Modes in MPM|thermodynamics modes]] for more details.
; <tt>double GetFnEnergy(double delta,double tr) const</tt><br><tt>double GetFtPrime(double delta,double tr) const</tt>
: These methods (one for normal and one for tangential directions) should return area under the traction law up to given displacement discontinuity in <tt>delta</tt> (in [[ConsistentUnits Command#Legacy and Consistent Units|Length units]]) at the calculated traction force in <tt>tr</tt> (in [[ConsistentUnits Command#Legacy and Consistent Units|Force units]]). The units of the return value should be [[ConsistentUnits Command#Legacy and Consistent Units|Energy units]]. Note that this current implementation of interfacial energy is based on assumption of an elastic imperfect interface. It can archived  in the [[MPM Global Archiving Options|global archiving interface energy option]], but has no other affect on calculations. It does not contribute to dissipated energy and does not cause any temperature changes. Future versions of imperfect interface laws may allow for dissipated energy calculations.

Latest revision as of 12:14, 14 April 2018

This section explains how to write C++ code create a custom imperfect interface for use in NairnMPM along with available imperfect interface contact laws.

Getting Started

The first steps are to create source files for the new interface law class, to give the law a unique name, and to allow the main NairnMPM code to instantiate the contact law to be used on cracks or in multimaterial mode.

Class Source Code

Custom imperfect interface should start by duplicating NonlinearInterface.cpp and NonlinearInterface.hpp files (or the files for one of its subclasses). All custom lows should inherit from the NonlinearInterface class. The new contact law will be used along with other Contact Laws and therefore will need a unique name and ID (an integer). In the new contact law's header file, replace NONLINEARINTERFACELAW copies above with a defined constant representing the new contact law's ID (which by convention is in UPPERCASE) and replace the number in the new constant's definition with the new ID. For example:

#define MYINTERFACE 201

All documented materials use numbers below 100. To avoid conflicts, those working on custom materials should use large numbers (>100).

Editing Required in Core Code

Almost all coding will be done in the new contact law class files, but for that law to be recognized as an option in NairnMPM, two places in the core code have to be edited first. These should be the only changes needed outside the new contact law's class files.

  1. Include the new contact law's header file at the top of Common/Read_XML/MaterialController.cpp:
     #include "Materials/MyInterface.hpp"
    

    or the appropriate relative path to the new contact law's header file.

  2. In MaterialController::AddMaterial(int matID,char *matName), add a new case in the switch(matID) section to call the default constructor of the new contact law when matID matches the new contact law's constant defined above.
    case MYINTERFACE:
       newMaterial = new MyInterface(matName,matID);
       break;
    

    If needed you can define a custom constructor and use that in this code (this need is rare).

Compiling with make Command

If you want to be able to do command-line compile using the make command, the new contact law source files need to be added to the file NairnMPM/build/makefile. The steps are best done by comparing to the NonlinearInterface entries in that makefile. The main makeFile editing steps are:

  1. Define an alias for relative path from makefile to the new material class files (in the list d)
    MyMaterial = $(src)/Materials/MyInterface
    
  2. Add an object file to the list of "all compiled objects" in the form Myintefacel.o.
  3. Add the new material's header file to the header file's needed to compile MaterialController.hpp. Using the alias defined in step 1, the added text will be $(MyInterfacel).hpp.
  4. Compile the new contact law's source file with lines such as
    MyInterfacl.o : $(MyInterface).cpp $(dprefix) $(MyInterface).hpp $(MaterialBase).hpp $(MPMBase).hpp \
    			$(CommonException).hpp
    	$(CC) $(CFLAGS) $(headers) -include $(prefix) $(MyInterface).cpp
    

    The list of headers must include all headers explicitly (or implicitly) included in the new contact law's source code. The leading white space for lines 2 (and on) must use only tabs and only one tab for the final compilation line.

Writing the Contact Law Class Source Code

The remaining steps can usually all be done by editing only the contact law's source code file. The requirements and optional methods are described in the following sections. This stage usually starts by keeping by deleting all methods in the copied file that are not listed below and then changing the class name in each definition (e.g., change NonlinearInterface:: to MyInterface:: in each method). Any methods deleted from the cpp file should have its template deleted from the hpp file as well.

Creating Material Properties

Usually a newly-created contact law will have additional material properties than the parent law's (although reuse parent law properties when possible). To create such properties, you need to define them, allow them in input command files (i.e., update the DTD file), set them in the contact law's class file, and use them in calculations. The following steps are needed:

  1. Define a class variable for the property in the .hpp file (usually int or double). It is best to define these properties in the protected section of the class, although public properties are sometimes needed (and are allowed).
  2. To allow the property in input files, select a unique property name (the name may or may not be same as variable in previous step). Define that property by name in the DTD file, usually as a simple XML element such as:
    <!ELEMENT prop (#PCDATA)>
    

    where prop is the new property name.

  3. Add that property's name to the list of allowed elements within the Material element definition (beginning in <!ELEMENT Material in the DTD file.
  4. Set some default value for the new property variable in the new material's constructor method.

One properties are signed, they are initialized and printed with the following methods:

Required Initialization Methods

MyInterface::MyInterface(char *matName,int matID) : NonlinearInterface(matName,matID)
This method is the default constructor. It can initialize any custom interface properties to default values. It should als call constructor for the parent contact law class. If no defaults are needed, this method can be omitted. The MyInterface::: is the custom class name. This term is omitted from subsequent methods names (but is needed).
char *InputContactProperty(char *xName,int &input,double &gScaling)
If xName string matches a property name for this material, set input to DOUBLE_NUM or INT_NUM (depending on the type of variable) and return a pointer to the class variable for that property. If gScaling is set, the input value will be multiplied by gScaling after it is read. If xName is not a recognized property return NonlinearInterface::InputMat(xName,input) to allow the super class to look for their valid property types (replace <tt NonlinearInterface with different name if you subclassed some other contact law).
void PrintContactLaw(void) const
This method should print all properties or call a super class and then print just the new law properties. Use a style similar to other contact laws. To help in formatting, you can use the MaterialBase class method PrintProperty(label,prop,units) to print a label, a property numeric value, and units within one column. You can use several calls in sequence to get several properties on the same line. You can also call PrintProperty(text,right) to print text in a column and right or left justified if right is true or false.

Required Accessors

const char *MaterialType(void)
Return a short name to describe the imperfect interface law. This string is printed with contact law properties in the output of simulations.

MPM Step Calculations

The remaining source code is to use the law in MPM calculations. Much of the work is done in the parent NonlinearInterface class. Your custom subclass only needs to return details about the contact law it self. The current implementation assumes a decoupled law that the normal and tangential forces are given by:

      [math]\displaystyle{ T_n = f_n([u_n]) \qquad {\rm and } \qquad T_t = f_t([u_t]) }[/math]

or where normal traction [math]\displaystyle{ T_n }[/math] is a function only of normal interfacial discontinuity [math]\displaystyle{ [u_n] }[/math] and tangential traction [math]\displaystyle{ T_t }[/math] is a function only of tangential interfacial discontinuity [math]\displaystyle{ [u_t] }[/math]. This limitation may change in the future. As a result, each method defined below has two version - one for normal tractions and one for tangential tractions.

Check for Stability

int CheckDnStability(double d,double m,double delta) const
int CheckDtStability(double d,double m,double delta) const
The size

Interface Law Properties

double GetFn(double delta,double area) const
double GetFt(double delta,double area) const
These methods (one for normal and one for tangential directions) should return the force due to the traction law at the provided the displacement discontinuity delta. The force must be in Force units for given displacement in Length units. Many imperfect interface laws give a traction or stress. To convert such a law to force, multiply the traction by the contact area (provided in the area).
double GetFnPrime(double delta,double area) const
double GetFtPrime(double delta,double area) const
These methods (one for normal and one for tangential directions) should return derivative of the force due to the traction law at the provided the displacement discontinuity delta. The derivative must be in Force units/Length units for given displacement in Length units. Many imperfect interface laws give a traction or stress which have derivate for stress per unit length. To convert such a law proper derivative units, multiply the traction derivative by the contact area (provided in the area).
double GetFnEnergy(double delta,double tr) const
double GetFtPrime(double delta,double tr) const
These methods (one for normal and one for tangential directions) should return area under the traction law up to given displacement discontinuity in delta (in Length units) at the calculated traction force in tr (in Force units). The units of the return value should be Energy units. Note that this current implementation of interfacial energy is based on assumption of an elastic imperfect interface. It can archived in the global archiving interface energy option, but has no other affect on calculations. It does not contribute to dissipated energy and does not cause any temperature changes. Future versions of imperfect interface laws may allow for dissipated energy calculations.