*****************************************************************************
                            Mike 1.0 README
                            4. March 2003
                    Discrete Element Code for Elle 
                    
                Copyright (c) D. Koehn and J. Arnold, 2003
*****************************************************************************

            Written by Daniel Koehn and Jochen Arnold
                    Oslo-Mainz 2002 to 2003
                big thanks to Anders Malthe-Srenssen
                    
                    koehn@mail.uni-mainz.de
                http://www.microstructure.info

This code is freeware and comes as it is. 
The authors make no representations about the suitability of this software for any purpose. 
The code contains bugs and is not written to be bullet proof. Please submit bug-reports 
to the above email address or to the Elle support site. 

All changes to the base code (except for the MyApplication.elle.cc files) should be recorded
in the version control file. Please submit changes to myself or to the elle site. See 
version control file for recent changes and fixed bugs. 
Please  comment any code that you add in detail ! 
(And remember to have fun !)


        Mainz, 4. March 2003 
        Daniel 
*******************************************************************************
                        INSTALLATION NOTICE
*******************************************************************************

    Installation notice: 
    change file unodes.h and unodes.cc in elle/basecode
    unodes.h -> add 
                        #include "unodesP.h"
    and add 
                    Unode *ElleGetParticleUnode(int id); 
                    
    to the list of functions. 
    unodes.cc -> add
    
    /*****************************************************************
    * Function that gets the id of a Unode in the vector *Unodes 
    * and returns a pointer of type Unode * to that Unode. 
    * Used by Mike.
    * Needed because cannot call to Unodes out of unodes.cc
    *
    * included Feb. 2002 daniel
    ****************************************************************/

    Unode *ElleGetParticleUnode(int id)
    {
    if (Unodes)  // if there is a vector of Unodes 
        {
        return &(*Unodes)[id]; // return pointer to that Unode
        }
    }
    
    to the file. Then you might have to reinstall. 
    
    probably will add this to the next official Elle distribution. 
    
    At the moment the imake file is not installed but we will probably do this
    in the version that can be downloaded. 
    Put Mike folder in the elle/processes folder. 
    Copy a compiled MakeFile into mike folder (for example gbm makefile)
    Go to line 100 and change: 
    
    CURRENT_DIR = processes/mike

    
    Go to line 1470 and change:
    
    SRCS=  mike.main.cc myApplication.elle.cc lattice.cc particle.cc 
    OBJS=  mike.main.o myApplication.elle.o lattice.o particle.o 
    THISAPP= myApplication
    
    now type make (you may have to run make depend first)
    
    and it should compile

    

    
******************************************************************************
                            HOW TO RUN MIKE
******************************************************************************

    INPUT:  to start mike run ./myApplication -i myInput.ppm.elle -f 1
    
    NOTE THAT MIKE NEEDS A VALID INPUT FILE GENERATED WITH THE NEW VERSION OF 
    PPM2ELLE. OTHERWISE IT WILL CRASH RIGHT AWAY !
    
    Switch out the unode triangulation, otherwise mike will crash when you try 
    to run it. 
    
    ________________________________________________________________________
    
    mike 1.0 has four main elements at the moment. 
    It is written in C++ mostly but contains C segments. 
    highest Class at the moment is the lattice class that builds everything and
    contains all the user function. Second class is the particle class that contains 
    single elastic particles. Lattice contains lists of particles, lists are circular 
    and can be run both ways. 
    Important for the USER are the user functions in the lattice class that are introduced
    below in detail. 
    the mike.main.cc file just defines some elle run time routines and is not so important. 
    The myApplication.elle.cc file contains the user stuff. This can be used like a normal 
    script where you just call your user function from the lattice class in order to do things
    in small loops etc.. 
    If you want to know more about the basic code go to the sources and look at the comments, 
    I tried to comment as much as I could. 
    
    ------------------------------------------------------------------------
    
    the myApplication.elle.cc file: 
    
    after the headers:
    
    Lattice *mike;  
    
    Defines a pointer of type object lattice called mike. You have to specify the name of your 
    object here (can call it anything you want). 
    
    in InitSetMike()     (this is the elle call function for the initialization of everything)
    
    mike = new Lattice;   
    
    This gives you the name of the object again and calls the constructor in the class 
    (builds basic lattice and connects to elle nodes, Grains(Flynns) etc. 
    
    After that you can call the initialization function in mike by calling: 
    
    mike->myInitializationFunction(myVariable);

    Functions in mike_1.0 are: 
    
    //------------------------------------------------------------------
    // functions for external use called from myApplication.elle.cc
    //------------------------------------------------------------------

    void WeakenGrain(int nb,float constant, float break_strength);   // changes grain constant and tensile strength

    void MakeGrainBoundaries(float constant, float break_strength); // changes grainboundary constant and tensile strength 

    void AdjustConstantGrainBoundaries();      // function to adjust modulus of grainboundary springs 

    void SetFracturePlot(int numbfrac, int local);     // plot after certain amount of fractures or not
 
    void ChangeParticle(int number, float shrink, int plot); // shrink or grow one particle
    
    void ShrinkGrain(int nb, float shrink, int plot);  // shrink or expand a whole grain

    void WeakenHorizontalLayer(double y_min,double y_max,float constant, float break_strength); // define layer 
 
    
    // sets a Gaussian distribution on spring constants of whole grains

    void SetGaussianSpringDistribution(double g_mean,double g_sigma);
 
    
    // sets pseudorandom distribution of spring constants of grains and of breaking strength of all springs

    void SetPhase(float str_mean, float str_size,float br_mean, float br_size);
    
    ------------------------------------------------------------------------
                            DETAILED DESCRIPTION
    ------------------------------------------------------------------------

    void WeakenGrain(int nb,float constant, float break_strength);   
    
        This function receives a number for a grain (integer), a number to change the spring constant of that 
        grain (float) and a number to change the breaking strength of the grain springs (float). Input numbers
        for constant or break_strength of 1.0 will change nothing i.e. the default value of the lattice will be 
        multiplied by this number. Therefore 2.0 will make springs 2 times as hard and 0.1 will make them a magnitude
        softer. 

    void MakeGrainBoundaries(float constant, float break_strength); 
    
        This function just receives a number to change the spring constant of grain Boundaries and a number to 
        change the breaking strength of grain boundaries. Also default values multiplied by this number (i.e. 1.0
        leaves default values). Can be used to lower breaking strength of all grain boundaries. 
    
     void AdjustConstantGrainBoundaries();    
     
        This function can be called to clean up grain strength settings. Normally the grain boundaries are not set 
        when you call WeakenGrain. Calling this function after you set some grain strength distribution will give
        grain boundary springs mean values of neighbouring grains. 

    void SetFracturePlot(int numbfrac, int local);    
    
        This little thing can be called at the beginning after the lattice 
        was made (in myApplication.elle.cc). First number is an integer after how many
        broken bonds the program will make a picture within deformation steps.  Is 
        useful sometimes when a lot of fractures develop after a deformation
        step and you want to see whats happening. Or you want to watch the 
        fracture develop. The second number is a number for the program to tell it
        that it should do a plot (put in 1 if you want pictures). 
        Default is 0 so if you dont call this function the 
        program will only make plots after time steps. 
 
    void ChangeParticle(int number, float shrink, int plot); 
    
        Change the size of particle number by amount shrink (can also expand depending on number, 
        radius of particle is mulitplied by this number). Variable plot = 0 is no picture and 1 
        is make a picture after the relaxation. Can also be called in the Run program.
    
    void ShrinkGrain(int nb, float shrink, int plot);  
  
        Change the size of a whole grain (int nb) by and amount shrink (again multiply by this number) and 
        plot a picture if variable plot is 1 and not if its 0. Can also be called in the Run program.
  
    void WeakenHorizontalLayer(double y_min,double y_max,float constant, float break_strength); 
 
        This function will change the spring constant or breaking strength of all grains that lie 
        within the boundaries y_min and y_max. This will make a horizontal layer with softer or harder
        grains.  Can be called several times. 
 

    void SetGaussianSpringDistribution(double g_mean,double g_sigma);
    
        This function sets a gaussian distribution for spring constants of grains. Input is the mean (center
        of distribution, absolute number) and the width of the distribution by sigma. Distribution is from 
        mean minus 2 sigma to mean plus 2 sigma. Default value for the spring constant in the model is 1.0. 
        This function generates a gaussian distribution of springconstants k_gauss for a 
        a given grain population. This means that the number of grains with a certain
        function sprinconstant f(k) within a desired intervall of values follows a  
        density function of the form f(k)=(sqrt(2*pi)*g_sigma)**(-1) exp((k_gauss-g_mean)/g_sigma)**2 
        where g_sigma is the standard variation, g_mean is the mean value.
        Remember that g_mean determines the centre of the distribution and 
        sigma the sharpness and height of the distribution. 
  
        the function receives the mean and the sigma for the distribution. 
        AS mentioned above the mean will be the mean of the spring constants of the 
        grains so that a mean statistic of the whole box should give this constant as 
        long as the system behaves linear elastic. The distibution is set for grains
        so that all particles in a grain will have the same spring constants (i.e. their
        springs will). mean will give absolute spring constants so that this function
        should be called first before anything else is set, otherwise it will be wiped out. 
        No problem however in setting a gauss distribution and afterwards making a layer
        hard or soft. Default spring constant in the model is 1.0 so that a mean of 1.0 
        will produce values clustering around 1.0 and 2.0 values with a mean spring constant 
        of 2.0 etc.. 
        We determine distribution within the interval mean plus minus two times sigma ! 
 
        The distribution is determined as follows: Pic a pseudorandom number (from mean 
        plus/minus 2 times sigma)  to give a grain
        a spring constant. Then determine the probability for that constant to appear from
        the Gauss function. Then pic a second pseudo random number (from 0 to sigma) to 
        determine whether or not the grain will actually get this spring constant. If 
        the random number is below the probability the spring contant is accepted, 
        if it is higher it is rejected and a new one is drawn. 
 
        Function writes also a text file with the set spring constant distribution
        calld gauss.txt

 
 
    void SetPhase(float str_mean, float str_size,float br_mean, float br_size);
    
        This function can be used to set a pseudoramdon distribution for spring constants of grains and for 
        all breaking strengths of all springs in the model. Function that sets a pseudorandom distribution on spring 	constants and/or on breaking-
        strength of springs. Spring constants are set for grains, breaking strength for the 
        whole model. First number is for the mean spring constant. This number is added or 
        substracted from the set values (default 1.0) so that 0.0 does not change the mean. 
        Therefore anisotropies or distributions in the model will not be distroyed but redistributed
        depending on what they are. The second number gives the distribution size. The distribution 
        is original mean times str_mean plus/minus half of str_size. 
        The same applies for the breaking strength of springs. 
        That means if the numbers are 0.0 0.4 0.8 0.6 and the mean spring constant was 1.2 and the 
        mean breaking strength 1.0 the spring constants are number from 1.0 to 1.4 and the 
        breaking strengths numbers from 1.7 to 2.3
        if second or fourth numbers are 0.0 no distributions are set for these values
        
        
    ----------------------------------------------------------------------------------------------
    
    SetMike() in the myApplication.elle.cc file is the runfunction for mike in elle. 
    
    You can call runfunction of mike in here in different loops etc. 
    
    ----------------------------------------------------------------------------------------------
                                    RUNFUNCTIONS
    ----------------------------------------------------------------------------------------------
    
    float DeformLattice(float move, int plot);     // uniaxial compression with side walls

    float DeformLatticePureShear(float move, int plot);   // deform volume constant pure shear

    void ChangeParticle(int number, float shrink, int plot); // shrink or grow one particle

    void ChangeParticleStress(float shrink, int plot); // shrink particle that is stressed most

    void ShrinkGrain(int nb, float shrink, int plot);  // shrink or expand a whole grain
    
    // function to dump statistics of a single grain
 
    void DumpStatisticStressGrain(double strain,int nb);
 
    // function to dump statistics of two grains

    void DumpStatisticStressTwoGrains(double strain,int nb1,int nb2);
 
    // function to dump statistics of a predifined box

    void DumpStatisticStressBox(double y_box_min,double y_box_max, double x_box_min,double x_box_max, double strain);
 
    -----------------------------------------------------------------------------------------------
                                DETAIL DESCRIPTION
    -----------------------------------------------------------------------------------------------
    Deformation: 
    -------------
    
    float DeformLattice(float move, int plot);  
    
        This function deforms the lattice by uniaxial movement in the y direction. The x walls are
        kept stationary. The upper wall is moved by a strain "move" in percent, positive is compression
        and negative tension. 1.0 would be 100 percent. Should be small value in the order of 0.0005. 
        int plot = 0 means no picture after relaxation and 1 means picture. The movements of particles in the
        box are averaged after each deformation step assuming that the particles are homogeneous. From this 
        condition the relaxation starts. This will produce no initial gradients in y in the box. 
        
    float DeformLatticePureShear(float move, int plot);
    
        This function deforms the lattice by pure shear area conservative deformation. move is again in percent, 
        where 1.0 is 100. plot is again 0 no pict and 1 make a pict. Moves upper wall downwards if move is positive
        and right wall to the right. Movements of particles are again averaged after a deformation step and then 
        relaxed so that no gradients develop initially. 
        
    Reaction/VolumeChange: 
    ----------------------
    
    void ChangeParticle(int number, float shrink, int plot); 
     
        This function changes the size of a particle (particle number by amount shrink). Plot again 0 no plot, 
        1 make plot. This way you can blow up particles.....  :-) :-) :-) 
        
    void ShrinkGrain(int nb, float shrink, int plot); 
    
        Yeah, this blows up whole grains ! Give it a grain number and a number to multiply the grain size and 
        the whole grain will do something (number larger 1.0 -> blow up, smaller -> shrink). plot again is 
        0 -> no plot and 1 -> plot picture. The number to change the volume of grains should be small because
        no average is taken so that this is quite something dramatic.... 

    void ChangeParticleStress(float shrink, int plot); 
    
        And this now includes the stress of the particle. I think its pressure, sxx + syy ( similar to the 
        mean stress) at the moment. It shrinks the particle with the highest stress by the amount shrink. 
        again plot 0 is no pict and plot 1 a pict. 
        This is not yet very advanced but can be changed and expanded easily. This is so far just to play 
        around. But can write lots of more functions based on the idea..... 
        
    Statistics:
    ------------
    
    void DumpStatisticStressGrain(double strain,int nb);
 
        Function that is used to dump a statistic file. This file gets the strain as 
        input from the main function (the strain that you put in the deformation file)
        and an int for the grain that you want the statistic for. Dumps mean values for 
        stress tensor and mean and differential stress for this grain. Called normally 
        by the user from the mike.elle.cc file within the Run section after deformation. 

        file is called statisticgrain.txt

        dumps strain, mean stress smax + smin/2, sxx, syy, sxy and differential stress
        smax - smin. Note that all compressive stresses are by definition negative ! 


    void DumpStatisticStressTwoGrains(double strain,int nb1,int nb2);
 
        ok, this is almost the same as the above code (to dump statistic for just
        one grain). This is now for two grains. So the function gets the strain 
        from the deformation you set and two grain numbers. This is useful to
        compare stresses in a hard layer and a soft layer. The function is an
        append file again, so the program writes the data directly to the file 
        each time and the file is there even if the code crashes (thanks, Jens)
        dumps number of grain, 
        strain on box, mean stress, differential stress and sxx, syy and sxy

    void DumpStatisticStressBox(double y_box_min,double y_box_max, double x_box_min,double x_box_max, double strain);

        Another file to dump statistics in a text file. This is for a whole box. It gets the input values
        box position with is box min in y box max in y, box min in x and box max in x and also the strain
        applied on the whole box (the elle box now) from your deformation file. 
        The function takes a mean of all particles in the box. It calculates means for the mean stress, 
        the differential stress and the stress tensor sxx, syy, sxy. 
        output is finite strain on box, mean stress, sxx,syy,sxy and differential stress. Function
        can be called from whereever you want it. 

                                


******************************************************************************

    About Mike: 
    
    Mike is a linear elastic discrete element code. 
    It contains circular particles that are connected by linear springs. 
    Springs can break once they exceed a breaking threshold (a tensile stress). 
    
    THE SCALE:
    
    Mike 1.0 uses a triangular lattice -> each particle is connected with six springs to neighbours
    (except for the boundaries). Note that this will produce anisotropies especially for fractures. 
    Distributions of breaking strength help to get rid of these directions. 
    
    Spring constants are by default 1.0. Breaking strength 0.0125. (set in particle.cc SetSprings()). 
    Youngs modulus is sqrt(3)/2 times spring constant. 
    
    The rest you scale yourself. If 1.0 in the model is 10 GPa for example then the breaking strength of 
    single grains is by default 125 MPa (you can change that by making grain boundaries weaker etc.), the youngs 
    modulus is sqrt(3)/2 * 10 GPa, any stress value in the statistic files is also value times 10000 in MPa.... 
    and so on. 
    
    The triangular lattice in 2d reproduces linear elasticity. 
    See analytical proof in:
    Flekky ...
    
    and additional literature in: 
    
    
    

    
