Creating MPM Materials

This section explains how to write C++ code create a new material class for use in NairnMPM.

Getting Started

The first steps are to create source files for the new material class, to give the material a unique name, and to allow the main NairnMPM code to instantiate the material class when it is needed for a particle.

Class Source Code

The best way to create the source code is to duplicate the NewMaterial.cpp and NewMaterial.hpp files, which are templates for a new material class. Edit the files and change NewMaterial to the new material's class name (e.g., MyMaterial). The NewMaterial template is a subclass of MaterialBase. Normally the new material will be a subclass of another class. If so, change MaterialBase references, as needed, to the actual parent class. These changes are needed whenever a method in the new material class needs to pass control to its immediate super class.

Material Class Hierarchy

The new material should be inserted into the NairnMPM material class hierarchy with its name (from previous step) and a unique ID (an integer). In the new material's header files, replace NEWMATERIAL with a defined constant representing the new material's ID (which by convention is in UPPERCASE) and replace the number in the new constant's definition with the new material's ID. For example:

#define MYMATERIAL 102

All documented materials use low numbers (starting with 1). 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 material class files, but for that material 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 material's class files.

Include the new material's header file at the top of Common/Read_XML/MaterialController.cpp:
```
 #include "Materials/MyMaterial.hpp"
```
or the appropriate relative path to the new material's header file.
In MaterialController::AddMaterial(int matID,char *matName), add a new case in the switch(matID) section to call the default constructor of the new material when matID matches the new material's constant defined above.
```
case MYMATERIAL:
   newMaterial=new MyMaterial(matName);
   break;
```
If needed you can define a custom constructor and use that in this code (this need is very rare).

Compiling with make Command

If you want to be able to do command-line compile using the make command, the new material source files needed to add to the file NairnMPM/build/makefile. The steps are best done by comparing to a similar material's files in that makefile. The main makeFile editing stepsare:

Define an alias for relative path from makefile to the new material class files (in the list d)
```
MyMaterial = $(src)/Materials/MyMaterial
```
Add an object file to the list of "all compiled objects" in the form MyMaterial.o.
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 $(MyMaterial).hpp.
Compile the new material's source file with lines such as
```
MyMaterial.o : $(MyMaterial).cpp $(dprefix) $(MyMaterial).hpp $(MaterialBase).hpp $(MPMBase).hpp \
			$(CommonException).hpp
	$(CC) $(CFLAGS) $(headers) -include $(prefix) $(MyMaterial).cpp
```
The list of headers must include all headers explicitly (or implicitly included in the new material'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 Material 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 virtual in the header file to insure the correct code is used at run time.

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 DTD file), set them in the material class file, validate them, and finally use them in calculations. The following steps are needed:

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).
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.
Add that property's name to the list of allowed elements within the Material element definition (beginning in <!ELEMENT Material in the DTD file.
Set some default value for the new property variable in the new material's constructor method.

Once a property variable and name are defined, you set that variable, check it setting, and use it with the following methods:

char *InputMat(char *xName,int &input): If xName string matches anew 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 xName is not a recognized property use return Parent::InputMat(xName,input) to allow the super class to look for their valid property types (replace Parent with the name of the superclass).

const char *VerifyAndLoadProperties(int np) (optional): 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 Parent::VerifyAndLoadProperties(np) to let super class check properties as well (replace Parent with the name of the superclass) The input np 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 ValidateForUse() 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).

FillTransportProperties(TransportProperties *) (optional): This method is called by MaterialBase::VerifyAndLoadProperties() and can be used calculate transport properties that are independent of the particle state and thus will remain constant throughout the calculation. The MaterialBase class automatically finds diffusion and conductivity tensors for isotropic materials by using material properties D (in diffusionCon) for diffusion constant and kCond (in kCond) for thermal conductivity. It only needs to be overridden is something more is needed.

void PrintMechanicalProperties(void): 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 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. This method need not pass on to MaterialBase because it does not print any mechanical properties. This method is called after VerifyAndLoadProperties() is done.

void PrintTransportProperties(void) (optional): This method should print all transport properties in format similar to other materials (see help PrintProperty() methods described above). It should only print them if transport is activated (i.e., if(DiffusionTask::active) or if(ConductionTask::active)). The MaterialBase class prints isotropic properties and the ORTHO and TRANSISO1(2) classes print anisotropic properties. No additional printing is needed if one of these classes handles the task.

void ValidateForUse(int np) (optional): This method is called just before the first MPM time step and only for materials used by one or more materials. Throw a CommonException if this material cannot by used in current MPM mode (type specified in np which will be PLANE_STRAIN_MPM, PLANE_STRESS_MPM, AXISYMMETRIC_MPM, or THREED_MPM). If no exceptions, must call the same method in the parent class. Basic material properties and usually checked in VerifyAndLoadProperties() instead; this method is for materials that may have valid properties, but may be contingent on other MPM settings.

Basic Class Accessors

These methods (required unless specified as optional) return basic facts about the material:

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

int MaterialTag(): Return the constant defined above.

double WaveSpeed(bool,MPMBase *): Return the maximum wave speed for material in mm/sec. 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.

double CurrentWaveSpeed(bool,MPMBase *) (optiona): Return the current wave speed, which might depend on particle state. This method is used by the AdjustTimeStep task to change time step as needed; without this method, the time step will not adjust for the material.

double ShearWaveSpeed(bool) (optional): Return the shear wave speed for the material. This method is only used by silent boundary conditions. The MaterialBase class returns WaveSpeed()/sqrt(3). You only need to override this method if you want a better value. Also note that silent boundary conditions only work for isotropic materials.

double MaximumDiffusion(void) (optional): Return maximum diffusion constant in cm²/sec. This method is called once for each material at the start of the calculations and after material properties are defined. The MaterialBase class returns the appropriate result for isotropic materials. You only need to override if you need a different result.

double MaximumDiffusivity(void) (optional): Return maximum thermal diffusivity in cm²/sec = 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 MaterialBase class returns the appropriate result for isotropic materials. You only need to override if you need a different result.

double GetCurrentRelativeVolume(MPMBase *) (optional): 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).

Tensor GetStress(Tensor *sp,double pressure) (optional): Return current stress. Most materials just return the contents of sp 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.

bool PartitionsElasticAndPlasticStrain(void) (optional): Materials that separate elastic and plastic strain should override and return TRUE. Otherwise the strain tensor on the particle is assumed to have the total strain and plastic strain can be used for anything else (e.g., hyperelastic materials track the elastic Left-Cauchy Green strain in the particle's plastic strain).

bool SupportsArtificialViscosity(void) (optional): Materials that implement artificial viscosity must override this method to return TRUE and then implement the calculations in the constitutive law methods.

History Dependent Properties

This section needs editing and include set initial properties method.

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:

char *MaterialData(void) (optional) - this method is called once for each material point using this material at the start of the calculations. It should allocate memory to store history data on the particle, initialize the values, and return a pointer to the structure. It can be as simple as a single variable or a pointer to an array of any number of variables. If a class overrides this method in a super class, the material will need to make sure all history data is correctly allocated and accessible. double GetHistory(int num,char *ptr) - this method should extract history variable number num from the allocated data pointer in ptr. This method is only used when archiving results and num will only be 1 through 4 (currently only the first four history variables can be archived). If num is out of range for the current material, this method should return zero. To access values of history variables during constitutive law calculations, you can fetch them from the material point methods GetHistoryPtr(), GetHistoryDble(), and SetHistoryDble().

MPM Step Calculations

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.

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

int SizeOfMechanicalProperties(int &altBufferSize) const: The size

void *GetCopyOfMechanicalProps(MPMBase *mptr,int np,void *matBuffer,void *altBuffer) const: 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 VerifyAndLoadProperties() instead. Some materials also used in FEA make use of LoadMechProps() in this task, but that method is limited to 2D because it only inputs a single rotation about the z axis.

void GetTransportProps(MPMBase *,int,TransportProperties *) 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 AddForces(). 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 diffusionTensor and kCondTensor variables. It is automatically handled for subclasses of TRANSISO1(2) if the only effect is the current orientation. Transport properties that never change (i.e., independent of particle state) should be calculated in FillTransportProperties() instead.

double GetHeatCapacity(MPMBase *) (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 using units of J/(g-K). This method can be omitted if heat capacity is a constant (note that the heatCapacity variable in MaterialBase is internally stored in J/(g-K) during calculations, although it is input as J/(kg-K)).

double GetCpMinusCv(MPMBase *mptr) (optional): Return the difference between constant-strain and constant-stress 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 J/(g-K). This call is only used by conduction calculations because heat flow is based on C_p instead of on C_v. For an elastic material, this difference is given by -M.αT/ρ, where M is stress-temperature tensor and α is thermal expansion tensor. For an isotropic material, this reduces to 3Kα²T/ρ where K is bulk modulus and α is linear thermal expansion coefficient. For an ideal gas, it reduces to nR/ρ. 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), one of these methods is called to implement the constitutive law:

void MPMConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np): 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, strain energy, plastic energy, heat energy, and any material history dependent variables defined for the material. To support thermal and solvent expansion, include their effect on the constitutive law. The input u is the velocity gradient times the time step, which gives displacement gradient for this time step. See note below on particle temperature updates.

void MPMConstLaw(MPMBase *,double,double,double,double,double,int): This method applies 2D constitutive law. It is deprecated, but still works; first method above is preferred.

void MPMConstLaw(MPMBase *,double,double,double,double,double,double,double,double,double,double,int): This method applies 3D constitutive law. It is deprecated, but still works; first method above is preferred.

These updates should never change the particle temperature, because that would invalidate some thermodynamics calculations in NairnMPM. To cause a particle temperature change, call IncrementHeatEnergy() 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.