User login

Navigation

You are here

Comments

Subscribe to Comments feed
Comments
Updated: 2 hours 4 min ago

A Twitter web on fracture

Mon, 2020-08-31 11:02

In reply to keep reading-13

Thank you for pointing to this paper. It looks fascinating, and I will study it.

In preparing for my graduate course on fracture, I have tweeted several threads: https://twitter.com/zhigangsuo/status/1300110928474705921. Hope they are useful to you.

Also, you will reach many readers on Twitter.

Re: poling

Mon, 2020-08-31 09:46

In reply to poling

Hi Sanjay,

For polling, I do it mainly through zoom. You can set up the poll (for a given zoom session) before the class in the meeting details on zoom. You can also do it live on the spot. I sometimes do the latter when a question comes up spontaneously in class. As the host, you will be able to see the results (and the names of students if you don't want it to be anonymous, e.g., for grading purposes). The polling is multiple choice, I usually end up doing T/F or Y/N. I mainly use this feature to add some level of engagement in the course and also to see if students are generally understanding the content. There are other polling services (e.g., pollev.com) where you can have more open-ended answering, i.e., students can submit a free response. I know some other instructors use this for group-work related activities. These resources have limitations though, e.g., the number of responses that you can receive is limited if you use the free version (hence again the better utility for group work).

For breakout sessions, yes, it is mainly for students to work on problems together and then come back to class. The students can “ask for help” in their given groups; I can see this as the host and then pop in whenever students have questions. Likewise, I can “move around” the various rooms to check in with each group. With large courses, moving around the groups can be tricky but the “ask for help” feature is helpful in this regard. I also have design projects in my classes, which the students do in groups. I schedule brainstorming sessions related to these projects during class time a few times each semester. Breakout sessions are helpful in this regard as well.

Thank you for the references

Mon, 2020-08-31 04:01

In reply to Pressure and twist

Dear Dr. Kolesnikov,

Thank you for these reference which were unknown to me so far. It seems highly relevant to the our study. 

Regards,

Nir

A twitter thread on fatigue-resistant polymer networks

Sun, 2020-08-30 12:23

In reply to Journal Club for March 2019: Fatigue of hydrogels

I am making a crosslink between this thread and a Twitter thread.

Thanks!

Sat, 2020-08-29 22:20

In reply to Fatigue-resistant polymer networks

Dear Zhigang,

Sure. Thanks a lot.  It really reminds me of the great experience of taking your fracture mechanics course 3 or 4 years ago!

Best,

Shaoting

live lecture

Sat, 2020-08-29 15:44

In reply to Last spring, I taught an

Ravi, I agree that the live lecture is important.  The regular schedule seems to help the students quite a bit.  I do record mine for everyone and find that they come to lecture and then watch it again later for select bits of information.

poling

Sat, 2020-08-29 15:40

In reply to Re: Remote Mechanics Teaching. How to do it well?

@Matt can you say more about how you design your polls.  Also for the breakouts, do you have them work on a question and then come back and report to the class?

Fatigue-resistant polymer networks

Sat, 2020-08-29 12:31

In reply to Engineering Sciences 247: Fracture Mechanics

I have just tweeted a thread on fatigue-resistant polymer networks.

Fatigue-resistant polymer networks

Sat, 2020-08-29 12:29

In reply to Journal Club for July 2020: Fatigue-resistant hydrogels: Principles, Experiments, and Applications

In preparing the graduate course on fracture, I have just tweeted a thread on fatigue-resistant polymer networks.

Dear Morteza, 

Fri, 2020-08-28 08:52

In reply to Discussion of fracture paper #26 - Cracks and anisotropic materials

Dear Morteza, 

Thank you so much for the clarification. I fully agree that the solution near the crack tip is the interesting part for everyone interested in putting the math to work. This is the meaning of life isn't it. 

It is interesting what you say about boundary collocation method. For a general application one would chose the Bousinesq-Cerruti solution that has no convergence  limitations, I guess. With a 2a long crack I would have tried (a2-z2)-1/2 times a polynomial with free coefficients. The z would definitely be z=x±iy for an isotropic material but for your anisotropic material the z=x±μy could be a good starting point. What are your thoughts? Per

Toughness of composites

Thu, 2020-08-27 08:09

In reply to Engineering Sciences 247: Fracture Mechanics

Twitter thread on toughness of composites

Fractocohesive length

Thu, 2020-08-27 07:55

In reply to Engineering Sciences 247: Fracture Mechanics

I have started a Twitter thread on fractocohesive length. Please join me for conversation.

Re: Remote Mechanics Teaching. How to do it well?

Wed, 2020-08-26 12:48

In reply to Remote Mechanics Teaching. How to do it well?

Great points made, Matt!

Like most people (I think) I have resisted the need to alter the traditional lecture format and have stuck with the methods that worked well when I was a student. However, the current situation has served as a positive nudge for me and I have learnt that online tools can indeed improve my teaching, if used correctly. 

This semester I am teaching two graduate classes--Theory of Elasticity and Structural Dynamics. Both classes are primarily online, though I hope to include some experiment lab sessions in the Structural Dynamics class. For both classes, I use Camtasia to record my lectures, using Onenote and my tablet to write the notes, and upload them as SCORM Contents on Blackboard. Camtasia lets me add quizzes at various locations in the videos; this keeps the students engaged, lets me reinforce important concepts, and gives me feedback on which concepts remain unclear for some students. Weekly assignment are posted and recieved on Blackboard as well. Some feedback that I got from students last semester was that they miss having the "structure" of needing to be in class at specific times and the discipline that forces upon them. So, I post the new videos for each week before the scheduled class hours--I am sure there is more I can do on this front. 

I am incorporating more face-to-face interaction through required Zoom meetings with the entire class one day per week, during prescheduled class hours. I use this session to revise the important concepts covered in previous weeks videos, given them problems to solve by themselves and then go over the solutions with them (this part I intend to move in-person when climate permits), let them ask any questions they might have, and quickly provide a quick introduction to what this week's videos will cover. As far as getting students to ask questions goes, I have found the same thing as Matt--students are more comfortable sending me a direct message using the chat feature than asking in front of everyone. I have not had any issues with students chatting with each other. 

One feature which I was missing was peer interaction as well as some level of interactive experience students get in a traditional classroom. It is difficult to recreate a "class-spirit". For this, I am experimenting with MS Teams. I have created a Team for the classroom with various channels related to class material, extra readings and resources, and just general interaction related to the subject. I post all important announcements here and encourage students to contact me directly using the chat function rather than emailing me. The chat within Teams lets me answer their questions more promptly and also allows us to do a quick video call and clarify things for them if we need to. I also keep posting any interesting related material I can find online (for eg. the EML webinars! Markus's webinar last week timed perfectly with me just finishing vibration and waves in strings--saved me a lecture recording since I asked them to watch it as part of their assignment!) and also encourage students to post anything they find useful and help each other discover new and cool things about the topic. I know one can do this in Blackboard too but it just doesn't work well there. Teams has worked out great for me, including for my research group. I also use the Class Onenote feature on Teams to share my written notes with the class.

Cheating is an issue for everyone, I am sure. I have taken an approach similar to Matt. Instead of trying to replicate in-class exams, I now give them more conceptual quizzes that are open-ended and force them to think and apply concepts covered in class. The quizzes are open-notes, open-book..I rely on them to not cheat off each other. I am more existential on that aspect--if they are intent on cheating then there is only so much I can do to prevent it. I trust that they are hurting themselves in the long run if they keep that attitude. On my part, I try to minimize opportunities for cheating and try to ensure an even playing field for all students without resorting to tools that force them to give up any semblance of privacy.

My wish is that there was one online platform that I could use for everything; that would keep things a bit easier for me and the students. But no such magic bullet exists so far, so I am trying to use the best aspects of the tools that are available to me. It has worked out quite well so far and student feedback has been quite positive. I know the approach needs more work, but my hope is to learn from this experience and add in-person components to only make things better. In the future, I am most excited about being able to use in-person time for problem solving rather than lecturing--truly flipping the class was a distant dream for me; turns out that it is easier than I thought it would be. 

Thanks for initiating this discussion, Sanjay. I would love to hear about and learn from other people's perspectives and how they are approaching this. 

Re: Remote Mechanics Teaching. How to do it well?

Tue, 2020-08-25 12:35

In reply to Remote Mechanics Teaching. How to do it well?

I have found pros and cons of teaching online. 

Pros:

  • With recorded lectures, students can watch at their own pace: rewind, pause, etc.
  • Some interactive features actually help facilitate participation, if done correctly. I have found that polling and breakout sessions are good for this purpose (although breakout sessions can be difficult to manage in large courses).
  • Related to the above bullet, one interesting thing I have found is that students participate even more when provided with the "chat" (text) feature in Zoom. I think they are less intimidated to use this feature (as opposed to raising their hand and talking in class) and have grown up texting and being online their entire lives. As a warning, enabling this feature can slow progress in your class and can be a distraction in trying to monitor it simultaneously while lecturing. I recommend that you tell the students only to use the chat feature to ask the instructor questions (not to communicate with each other). Even though it slows the class down, to me it is worth it in terms of the benefits of increased participation in class.

Cons:

  • Last semester I uploaded pre-recorded lectures and then held online office hours for discusison sessions. I found that very few students attended office hours (whereas 90%+ came to class when it occurred in person).
  • Related to the above bullet, I think it is tempting for students to put off watching pre-recorded lectures (e.g., and try to watch a bunch of them all at once near the end of the semester). Likewise, I am sure they are watching tv, their phones, another copmuter monitor, etc. while "watching" the lectures. With their attention divided, they will not retain information as well.
  • Cheating. A whole other can of worms . . . I try to give open-ended, conceptual questions to help here.

 

Overall, if we have to teach online, I agree with Ravi (and his students) that doing the lecture live is still very important. It gives students a structured time to attend the course and ask questions. Incorporating some interactive components, even if they seem trivial (breakout sessions in zoom, polling in zoom (even Y/N or T/F), short quizzes during class) is vital, particularly in large classes, to keep the students' attention. An even better version would be to use some interactive "learning modules", e.g., through augmented reality and virtual reality. I have attempted to do this myself and through hiring some undergraduate/graduate students to help (from CS) but implementation is tricky. The know-how to develop these modules requires more time than most faculty members are willing to commit. I have also worked with a couple of AR/VR companies in this regard. Partnering with a company greatly decreases the time for implementation and increases the quality of the final product. The issue there comes in through user rights; typically the companies want to charge per user per semester, which is not sustainable for a single instructor and perhaps not even for a large department.

Last spring, I taught an

Mon, 2020-08-24 13:09

In reply to Remote Mechanics Teaching. How to do it well?

Last spring, I taught an undergrad mechanics of materials class and a graduate fracture mechanics class. The undergrad class had 160 students. Prior to transition online, I wrote things on an iPad using GoodNotes; the lectures were recorded and available to students for asynchronous viewing. So, the only thing that changed was that instead of in-person delivery, I recorded the lectures ahead of time and released one at each assigned lecture time. Students commented that while the transition worked well, the lack of interaction was an issue - they felt out of touch. The graudate class had 10 students; after the very first recorded lecture, they indicated that they prefer live lectures, and so that is what I did, still writing on the iPad. I would always preload into my iPad, all the files, figures, problem statements, etc necessary for each class so that all materials are accessed through a single device. 

This fall, I am scheduled to teach Introduction of Solid Mechanics (Linear); I plan to do live lectures with the iPad-zoom interface.

I think the overall feedback about the online delivery is the absence of connections. So, this fall I am planning to hold a town-hall once a week for 30 minutes or so and provide an opportunity for additional interaction with and between students.

Innovations? I am thinking of animations similar to gaming, but this is resource heavy!

 

Things get interesting in mechanotransduction

Sun, 2020-08-23 08:16

In reply to Journal Club for June 2020: Mechanically instructive biomaterials: a synergy of mechanics, materials and biology

Dear Zhenwei Ma and Jianyu Li, 

Thanks for taking up the discussion on adhesion and moving towards mechanotransduction. I really enjoyed your post. 

For a cell to probe the mechanical properties of the matrix, integrins first attach to the matrix ligand(adhesion). For example, integrins can bind to RGD of the matrix. Then they probe the matrix mechanical properties at a some frequency. Depending on the properties they sense, integrins can cluster and a cascade of signalling events occur inside the cell all the way up to translocation of YAP proteins inside the nucleus which can affect the genes, etc. Now with viscoelastic matrices, things get complicated. To make it simple, people used stiffness as a paramter and studied how cell behavior changes with stiffness. As you pointed out, Mooney and Chaudhuri's group worked on time dependent aspects. Also recently Paul Janmey's group works on viscoelastic solid like matrices. In addtion to this, there are papers which say the pore size of the gel determines cellular mechanotransduction. This is because the pores present adhesive ligands (RGD) at different  spacings. Also, the beauty of mechanotransduction is it is is biphasic. Matrix affects cell behavior -> also cells can degrade the matrix. For example, in Ovijit chaudhuri's group one of the reason they say for differences in cell spreading is that, cells remodel the matrix when you have a relaxing matrix. 

Now coming back to adhesion, doesnt adhesion would mean different things for cells and tissues, isnt it? Does a highly adhesive material would really matter to cells since they care more about the underlying mechanical properties and they anyway want to attach and spread to the matrix ligands? Or am I missing something here?

In this case, how do you think we can use tissue adhesives for cellular mechanotransduction studies? 

RE:VUMAT

Sat, 2020-08-22 19:05

In reply to Abaqus mailing list

I apologize for the inconvenience, but will it be possible to send the VUMAT to my email address? rctron@banetc.com

Abaqus mailing list

Sat, 2020-08-22 17:12

In reply to VUMAT implementation of crystal plasticity

subscribe to and seek assistance from the
ABAQUS mailing list
https://groups.yahoo.com/neo/groups/Abaqus/info
Search the archive of the list before posting in it.

RE:VUMAT

Sat, 2020-08-22 14:02

In reply to VUMAT

Thank you so much for the VUMAT.  However, I get "The executable package.exe aborted with system error code 1073740940" when I run the subroutine. What can be the reason for this error?

VUMAT

Sat, 2020-08-22 07:35

In reply to VUMAT implementation of crystal plasticity

Hello,

I obtained the code attached underneath from someone in a similar discussion. My browser does not show an option to upload a file. The source of the code was not disclosed to me.

Other possibility:
Interface between umat (Abaqus/Standard) and vumat (Abaqus/Explicit).
https://soilmodels.com/download/vumat_umatinterface-zip/

Either way: no guarantee from my side.

Good luck,

Frank

 

-------------------------------------------------------------------------------

 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C *** Implicit CPFE Vumat
C
C *** Modifications:
C
C    09/03/2009 by S Wang
C
C    1: Stress component order: 11, 22, 33, 12, 23, 31
C    2: Sequence of Slpdef, Slpspn
C    3: Ddemsd(direct component corrected)
C    4: Shear stran components stored in tensor component, not engineering compoents
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine vumat (
     1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2 stepTime, totalTime, dt, cmname, coordMp, charLength,
     3 props, density, strainInc, relSpinInc,
     4 tempOld, stretchOld, defgradOld, fieldOld,
     5 stressOld, stateOld, enerInternOld, enerInelasOld,
     6 tempNew, stretchNew, defgradNew, fieldNew,
     7 stressNew, stateNew, enerInternNew, enerInelasNew )

c-------------------------------------------------------------------------------
c  (1)  The list of variables above defining the *UMAT subroutine,
c  and the first (standard) block of variables dimensioned below,
c  have variable names added compared to earlier ABAQUS versions.
c
c  (2)  The statement: include 'aba_param.inc' must be added as below.
c
c  (3)  ABAQUS files use double precision only.
c    The file aba_param.inc has a line "implicit real*8" and, since
c    it is included in the main subroutine, it will define the variables
c    there as double precision.  But other subroutines still need the
c    definition "implicit real*8" since there may be variables that are
c    not passed to them through the list or common block.
c
c  (4)  This is for ABAQUS version 6 or newer
c
c  (5)  This UMAT has been implemented with the Bassani and Wu hardening law.
c  
c   The UMAT has been fixed by adding state variables to keep track of the
c   *total* slip on each slip system by integrating up the absolute value
c   of slip rates for each individual slip system.  These "solution dependent
c   variables" are available for postprocessing.  The only required change
c   in the input file is that the DEPVAR command must be changed.
c
C-----  Subroutines:
C
C       ROTATION     -- forming rotation matrix, i.e. the direction
C                       cosines of cubic crystal [100], [010] and [001]
C                       directions in global system at the initial
C                       state
C
C       SLIPSYS      -- calculating number of slip systems, unit
C                       vectors in slip directions and unit normals to
C                       slip planes in a cubic crystal at the initial
C                       state
C
C       GSLPINIT     -- calculating initial value of current strengths
C                       at initial state
C
C       STRAINRATE   -- based on current values of resolved shear
C                       stresses and current strength, calculating
C                       shear strain-rates in slip systems
C
C       LATENTHARDEN -- forming self- and latent-hardening matrix
C
C       ITERATION    -- generating arrays for the Newton-Rhapson
C                       iteration
C
C       LUDCMP       -- LU decomposition
C
C       LUBKSB       -- linear equation solver based on LU
C                       decomposition method (must call LUDCMP first)

c-------------------------------------------------------------------------------

C-----  Variables:
C
C       stressOld, stressNew,  -- stresses (INPUT & OUTPUT)
C       stateOld, stateNew, -- solution dependent state variables (INPUT & OUTPUT)

C-----  Variables passed in for information:
C
C       DSTRAN -- increments of strains
C       CMNAME -- name given in the *MATERIAL option
C       NDIR   -- number of direct stress components (sigma_xx,sigma_yy,sigma_zz)
C       NSHR   -- number of engineering shear stress components (sigma_yx,sigma_zx,sigma_yz)
C       ndir+nshr  -- NDIR+NSHR
C       nstatev -- number of solution dependent state variables (as
C                 defined in the *DEPVAR option)
C       PROPS  -- material constants entered in the *USER MATERIAL
C                 option
C       NPROPS -- number of material constants
C

C-----  This subroutine provides the plastic constitutive relation of
C     single crystals for finite element code ABAQUS.  The plastic slip
C     of single crystal obeys the Schmid law.  The program gives the
C     choice of small deformation theory and theory of finite rotation
C     and finite strain.
C       The strain increment is composed of elastic part and plastic
C     part.  The elastic strain increment corresponds to lattice
C     stretching, the plastic part is the sum over all slip systems of
C     plastic slip.  The shear strain increment for each slip system is
C     assumed a function of the ratio of corresponding resolved shear
C     stress over current strength, and of the time step.  The resolved
C     shear stress is the double product of stress tensor with the slip
C     deformation tensor (Schmid factor), and the increment of current
C     strength is related to shear strain increments over all slip
C     systems through self- and latent-hardening functions.

C-----  The implicit integration method proposed by Peirce, Shih and
C     Needleman (1984) is used here.  The subroutine provides an option
C     of iteration to solve stresses and solution dependent state
C     variables within each increment.

C-----  The present program is for a single CUBIC crystal.  However,
C     this code can be generalized for other crystals (e.g. HCP,
C     Tetragonal, Orthotropic, etc.).  Only subroutines ROTATION and
C     SLIPSYS need to be modified to include the effect of crystal
C     aspect ratio.
C

C-----  Important notice:
C
C     (1) The number of state variables nstatev must be larger than (or
CFIX      equal to) TEN (10) times the total number of slip systems in
C         all sets, NSLPTL, plus FIVE (5)
CFIX           nstatev >= 10 * NSLPTL + 5
C         Denote s as a slip direction and m as normal to a slip plane.
C         Here (s,-m), (-s,m) and (-s,-m) are NOT considered
C         independent of (s,m).  The number of slip systems in each set
C         could be either 6, 12, 24 or 48 for a cubic crystal, e.g. 12
C         for {110}<111>.
C
C         Users who need more parameters to characterize the
C         constitutive law of single crystal, e.g. the framework
C         proposed by Zarka, should make nstatev larger than (or equal
C         to) the number of those parameters NPARMT plus nine times
C         the total number of slip systems, NSLPTL, plus five
CFIX           nstatev >= NPARMT + 10 * NSLPTL + 5
C
C     (2) The tangent stiffness matrix in general is not symmetric if
C         latent hardening is considered.  Users must declare "UNSYMM"
C         in the input file, at the *USER MATERIAL card.
C
c-------------------------------------------------------------------------------
c
    INCLUDE 'VABA_PARAM.INC'

      CHARACTER*8 CMNAME
    Parameter    (ND=150)

C-----  The parameter ND determines the dimensions of the arrays in
C     this subroutine.  The current choice 150 is a upper bound for a
C     cubic crystal with up to three sets of slip systems activated. 
C     Users may reduce the parameter ND to any number as long as larger
C     than or equal to the total number of slip systems in all sets. 
C     For example, if {110}<111> is the only set of slip system
C     potentially activated, ND could be taken as twelve (12). 

C    FOR VUMAT:
C    All arrays dimensioned as (*) are not used in this subroutine
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     2 charLength(*), strainInc(nblock,ndir+nshr),
     3 relSpinInc(nblock,nshr+nshr), tempOld(*),
     4 stretchOld(*),defgradOld(*),
     5 fieldOld(*), stressOld(nblock,ndir+nshr),
     6 stateOld(nblock,nstatev), enerInternOld(nblock),
     7 enerInelasOld(nblock), tempNew(*),
     8 stretchNew(*),defgradNew(*),
     9 fieldNew(*),stressNew(nblock,ndir+nshr),
     + stateNew(nblock,nstatev),enerInternNew(nblock),
     + enerInelasNew(nblock)

C    FOR UMAT:
      dimension stress(ndir+nshr),statev(nstatev),dstran(ndir+nshr)
c
      DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),
     2          SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),
     3          SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),
     4          DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),
     5          FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),
     6          H(ND,ND), 
     7          DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),
     8          DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),
     9          WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)

      DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),
     2          GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),
     3          DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),
     4          DSPNRO(3,ND), DSPDRO(3,ND),
     5          DHDGDG(ND,ND)

C-----  NSLIP  -- number of slip systems in each set
C-----  SLPDIR -- slip directions (unit vectors in the initial state)
C-----  SLPNOR -- normals to slip planes (unit normals in the initial
C                 state)
C-----  SLPDEF -- slip deformation tensors (Schmid factors)
C                 SLPDEF(1,i) -- SLPDIR(1,i)*SLPNOR(1,i)
C                 SLPDEF(2,i) -- SLPDIR(2,i)*SLPNOR(2,i)
C                 SLPDEF(3,i) -- SLPDIR(3,i)*SLPNOR(3,i)
C                 SLPDEF(4,i) -- SLPDIR(1,i)*SLPNOR(2,i)+
C                                SLPDIR(2,i)*SLPNOR(1,i)
C                 SLPDEF(5,i) -- SLPDIR(2,i)*SLPNOR(3,i)+
C                                SLPDIR(3,i)*SLPNOR(2,i)
C                 SLPDEF(6,i) -- SLPDIR(1,i)*SLPNOR(3,i)+
C                                SLPDIR(3,i)*SLPNOR(1,i)

C                 where index i corresponds to the ith slip system
C-----  SLPSPN -- slip spin tensors (only needed for finite rotation)
C                 SLPSPN(1,i) -- [SLPDIR(1,i)*SLPNOR(2,i)-
C                                 SLPDIR(2,i)*SLPNOR(1,i)]/2
C                 SLPSPN(2,i) -- [SLPDIR(2,i)*SLPNOR(3,i)-
C                                 SLPDIR(3,i)*SLPNOR(2,i)]/2
C                 SLPSPN(3,i) -- [SLPDIR(3,i)*SLPNOR(1,i)-
C                                 SLPDIR(1,i)*SLPNOR(3,i)]/2

C                 where index i corresponds to the ith slip system
C-----  DSPDIR -- increments of slip directions
C-----  DSPNOR -- increments of normals to slip planes
C
C-----  DLOCAL -- elastic matrix in local cubic crystal system
C-----  D      -- elastic matrix in global system
C-----  ROTD   -- rotation matrix transforming DLOCAL to D
C
C-----  ROTATE -- rotation matrix, direction cosines of [100], [010]
C                 and [001] of cubic crystal in global system
C
C-----  FSLIP  -- shear strain-rates in slip systems
C-----  DFDXSP -- derivatives of FSLIP w.r.t x=TAUSLP/GSLIP, where
C                 TAUSLP is the resolved shear stress and GSLIP is the
C                 current strength
C
C-----  DDEMSD -- double dot product of the elastic moduli tensor with
C                 the slip deformation tensor plus, only for finite
C                 rotation, the dot product of slip spin tensor with
C                 the stress
C
C-----  H      -- self- and latent-hardening matrix
C                 H(i,i) -- self hardening modulus of the ith slip
C                           system (no sum over i)
C                 H(i,j) -- latent hardening molulus of the ith slip
C                           system due to a slip in the jth slip system
C                           (i not equal j)
C
C-----  DSTRES -- Jaumann increments of stresses, i.e. corotational
C                 stress-increments formed on axes spinning with the
C                 material
C-----  DELATS -- strain-increments associated with lattice stretching
C                 DELATS(1) - DELATS(3) -- normal strain increments
C                 DELATS(4) - DELATS(6) -- engineering shear strain
C                                          increments
C-----  DSPIN  -- spin-increments associated with the material element
C                 DSPIN(1) -- component 12 of the spin tensor
C                 DSPIN(2) -- component 31 of the spin tensor
C                 DSPIN(3) -- component 23 of the spin tensor
C
C-----  DVGRAD -- increments of deformation gradient in the current
C                 state, i.e. velocity gradient times the increment of
C                 time
C
C-----  DGAMMA -- increment of shear strains in slip systems
C-----  DTAUSP -- increment of resolved shear stresses in slip systems
C-----  DGSLIP -- increment of current strengths in slip systems
C
C
C-----  Arrays for iteration:
C
C            FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,
C            DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO,
C            DHDGDG
C
C
C-----  Solution dependent state variable STATEV:
C            Denote the number of total slip systems by NSLPTL, which
C            will be calculated in this code.
C
C       Array STATEV:
C       1          - NSLPTL    :  current strength in slip systems
C       NSLPTL+1   - 2*NSLPTL  :  shear strain in slip systems
C       2*NSLPTL+1 - 3*NSLPTL  :  resolved shear stress in slip systems
C
C       3*NSLPTL+1 - 6*NSLPTL  :  current components of normals to slip
C                                 slip planes
C       6*NSLPTL+1 - 9*NSLPTL  :  current components of slip directions
C
CFIX    9*NSLPTL+1 - 10*NSLPTL :  total cumulative shear strain on each
CFIX                              slip system (sum of the absolute
CFIX                              values of shear strains in each slip
CFIX                              system individually)
CFIX
CFIX    10*NSLPTL+1             : total cumulative shear strain on all
C                                 slip systems (sum of the absolute
C                                 values of shear strains in all slip
C                                 systems)
c
C #######    ---change ---nstatev from 125 to 135  #########
C
CFIX    10*NSLPTL+2 - nstatev-4-10 : additional parameters users may need
C                                  to characterize the constitutive law
C                                  of a single crystal (if there are
C                                  any).
C
C       nstatev-3-10               :  number of slip systems in the 1st set
C       nstatev-2-10               :  number of slip systems in the 2nd set
C       nstatev-1-10               :  number of slip systems in the 3rd set
C       nstatev  -10               :  total number of slip systems in all
C                                 sets
C
C
C-----  Material constants PROPS:
C
C       PROPS(1) - PROPS(21) -- elastic constants for a general elastic
C                               anisotropic material
C
C            isotropic   : PROPS(i)=0  for  i>2
C                          PROPS(1) -- Young's modulus
C                          PROPS(2) -- Poisson's ratio
C
C            cubic       : PROPS(i)=0  for i>3
C                          PROPS(1) -- c11
C                          PROPS(2) -- c12
C                          PROPS(3) -- c44
C
C            orthotropic : PORPS(i)=0  for  i>9
C                          PROPS(1) - PROPS(9) are D1111, D1122, D2222,
C                          D1133, D2233, D3333, D1212, D1313, D2323,
C                          respectively, which has the same definition
C                          as ABAQUS for orthotropic materials
C                          (see *ELASTIC card)
C
C            anisotropic : PROPS(1) - PROPS(21) are D1111, D1122,
C                          D2222, D1133, D2233, D3333, D1112, D2212,
C                          D3312, D1212, D1113, D2213, D3313, D1213,
C                          D1313, D1123, D2223, D3323, D1223, D1323,
C                          D2323, respectively, which has the same
C                          definition as ABAQUS for anisotropic
C                          materials (see *ELASTIC card)
C
C
C       PROPS(25) - PROPS(56) -- parameters characterizing all slip
C                                systems to be activated in a cubic
C                                crystal

C
C            PROPS(25) -- number of sets of slip systems (maximum 3),
C                         e.g. (110)[1-11] and (101)[11-1] are in the
C                         same set of slip systems, (110)[1-11] and
C                         (121)[1-11] belong to different sets of slip
C                         systems
C                         (It must be a real number, e.g. 3., not 3 !)
C
C            PROPS(33) - PROPS(35) -- normal to a typical slip plane in
C                                     the first set of slip systems,
C                                     e.g. (1 1 0)
C                                     (They must be real numbers, e.g.
C                                      1. 1. 0., not 1 1 0 !)
C            PROPS(36) - PROPS(38) -- a typical slip direction in the
C                                     first set of slip systems, e.g.
C                                     [1 1 1]
C                                     (They must be real numbers, e.g.
C                                      1. 1. 1., not 1 1 1 !)
C
C            PROPS(41) - PROPS(43) -- normal to a typical slip plane in
C                                     the second set of slip systems
C                                     (real numbers)
C            PROPS(44) - PROPS(46) -- a typical slip direction in the
C                                     second set of slip systems
C                                     (real numbers)
C
C            PROPS(49) - PROPS(51) -- normal to a typical slip plane in
C                                     the third set of slip systems
C                                     (real numbers)
C            PROPS(52) - PROPS(54) -- a typical slip direction in the
C                                     third set of slip systems
C                                     (real numbers)
C
C
C        PROPS(57) - PROPS(72) -- parameters characterizing the initial
C                                orientation of a single crystal in
C                                global system
C            The directions in global system and directions in local
C            cubic crystal system of two nonparallel vectors are needed
C            to determine the crystal orientation.
C
C            PROPS(57) - PROPS(59) -- [p1 p2 p3], direction of first
C                                     vector in local cubic crystal
C                                     system, e.g. [1 1 0]
C                                     (They must be real numbers, e.g.
C                                      1. 1. 0., not 1 1 0 !)
C            PROPS(60) - PROPS(62) -- [P1 P2 P3], direction of first
C                                     vector in global system, e.g.
C                                     [2. 1. 0.]
C                                     (It does not have to be a unit
C                                      vector)
C
C            PROPS(65) - PROPS(67) -- direction of second vector in
C                                     local cubic crystal system (real
C                                     numbers)
C            PROPS(68) - PROPS(70) -- direction of second vector in
C                                     global system
C
C
C       PROPS(73) - PROPS(96) -- parameters characterizing the visco-
C                                plastic constitutive law (shear
C                                strain-rate vs. resolved shear
C                                stress), e.g. a power-law relation
C
C            PROPS(73) - PROPS(80) -- parameters for the first set of
C                                     slip systems
C            PROPS(81) - PROPS(88) -- parameters for the second set of
C                                     slip systems
C            PROPS(89) - PROPS(96) -- parameters for the third set of
C                                     slip systems
C
C
C       PROPS(97) - PROPS(144)-- parameters characterizing the self-
C                                and latent-hardening laws of slip
C                                systems
C
C            PROPS(97) - PROPS(104)-- self-hardening parameters for the
C                                     first set of slip systems
C            PROPS(105)- PROPS(112)-- latent-hardening parameters for
C                                     the first set of slip systems and
C                                     interaction with other sets of
C                                     slip systems
C
C            PROPS(113)- PROPS(120)-- self-hardening parameters for the
C                                     second set of slip systems
C            PROPS(121)- PROPS(128)-- latent-hardening parameters for
C                                     the second set of slip systems
C                                     and interaction with other sets
C                                     of slip systems
C
C            PROPS(129)- PROPS(136)-- self-hardening parameters for the
C                                     third set of slip systems
C            PROPS(137)- PROPS(144)-- latent-hardening parameters for
C                                     the third set of slip systems and
C                                     interaction with other sets of
C                                     slip systems
C
C
C       PROPS(145)- PROPS(152)-- parameters characterizing forward time
C                                integration scheme and finite
C                                deformation
C
C            PROPS(145) -- parameter theta controlling the implicit
C                          integration, which is between 0 and 1
C                          0.  : explicit integration
C                          0.5 : recommended value
C                          1.  : fully implicit integration
C
C            PROPS(146) -- parameter NLGEOM controlling whether the
C                          effect of finite rotation and finite strain
C                          of crystal is considered,
C                          0.        : small deformation theory
C                          otherwise : theory of finite rotation and
C                                      finite strain
C
C
C       PROPS(153)- PROPS(160)-- parameters characterizing iteration
C                                method
C
C            PROPS(153) -- parameter ITRATN controlling whether the
C                          iteration method is used,
C                          0.        : no iteration
C                          otherwise : iteration
C
C            PROPS(154) -- maximum number of iteration ITRMAX
C
C            PROPS(155) -- absolute error of shear strains in slip
C                          systems GAMERR
C

    DO 100 NK9=1,NBLOCK
c
c
        DO NK8=1,ndir+nshr
        STRESS(NK8)=STRESSOLD(NK9,NK8)
        DSTRAN(NK8)=STRAININC(NK9,NK8)
      END DO

      DO NK8=1,nstatev
        STATEV(NK8)=STATEOLD(NK9,NK8)
      END DO
c
C-----  Elastic matrix in local cubic crystal system: DLOCAL
c
      DO J=1,6
         DO I=1,6
            DLOCAL(I,J)=0.
         END DO
      END DO

      CHECK=0.
      DO J=10,21
         CHECK=CHECK+ABS(PROPS(J))
      END DO

      IF (CHECK.EQ.0.) THEN
         DO J=4,9
            CHECK=CHECK+ABS(PROPS(J))
          END DO

         IF (CHECK.EQ.0.) THEN

            IF (PROPS(3).EQ.0.) THEN
c
C-----  Isotropic material
c
               GSHEAR=PROPS(1)/(1.+PROPS(2))
               E11=2.*GSHEAR*(1.-PROPS(2))/(1.-2.*PROPS(2))
               E12=2.*GSHEAR*PROPS(2)/(1.-2.*PROPS(2))

               DO J=1,3
                  DLOCAL(J,J)=E11

                  DO I=1,3
                     IF (I.NE.J) DLOCAL(I,J)=E12
                  END DO

                  DLOCAL(J+3,J+3)=GSHEAR
               END DO

            ELSE

C-----  Cubic material

               DO J=1,3
                  DLOCAL(J,J)=PROPS(1)

                  DO I=1,3
                     IF (I.NE.J) DLOCAL(I,J)=PROPS(2)
                  END DO

                  DLOCAL(J+3,J+3)=PROPS(3)
               END DO

            END IF

         ELSE

C-----  Orthotropic metarial

            DLOCAL(1,1)=PROPS(1)
            DLOCAL(1,2)=PROPS(2)
            DLOCAL(2,1)=PROPS(2)
            DLOCAL(2,2)=PROPS(3)

            DLOCAL(1,3)=PROPS(4)
            DLOCAL(3,1)=PROPS(4)
            DLOCAL(2,3)=PROPS(5)
            DLOCAL(3,2)=PROPS(5)
            DLOCAL(3,3)=PROPS(6)

            DLOCAL(4,4)=PROPS(7)
            DLOCAL(5,5)=PROPS(8)
            DLOCAL(6,6)=PROPS(9)

         END IF

      ELSE

C-----  General anisotropic material

         ID=0
         DO J=1,6
            DO I=1,J
               ID=ID+1
               DLOCAL(I,J)=PROPS(ID)
               DLOCAL(J,I)=DLOCAL(I,J)
            END DO
         END DO
      END IF

C-----Rotation matrix: ROTATE, i.e. direction cosines of [100], [010]
C     and [001] of a cubic crystal in global system
C
      CALL ROTATION (PROPS(57), ROTATE)

C-----  Rotation matrix: ROTD to transform local elastic matrix DLOCAL
C     to global elastic matrix D
C
      DO J=1,3
         J1=1+J/3
         J2=2+J/2

         DO I=1,3
            I1=1+I/3
            I2=2+I/2

            ROTD(I,J)=ROTATE(I,J)**2
            ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTATE(I,J2)
            ROTD(I+3,J)=ROTATE(I1,J)*ROTATE(I2,J)
            ROTD(I+3,J+3)=ROTATE(I1,J1)*ROTATE(I2,J2)+
     2                    ROTATE(I1,J2)*ROTATE(I2,J1)

         END DO
      END DO

C-----Elastic matrix in global system: D
C     {D} = {ROTD} * {DLOCAL} * {ROTD}transpose
C
      DO J=1,6
         DO I=1,6
            D(I,J)=0.
         END DO
      END DO

      DO J=1,6
         DO I=1,J

            DO K=1,6
               DO L=1,6
                  D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L)
               END DO
            END DO

            D(J,I)=D(I,J)

         END DO
      END DO

C-----  Total number of sets of slip systems: NSET
      NSET=NINT(PROPS(25))
      IF (NSET.LT.1) THEN
         WRITE (6,*) '***ERROR - zero sets of slip systems'
         STOP
      ELSE IF (NSET.GT.3) THEN
         WRITE (6,*)
     2     '***ERROR - more than three sets of slip systems'
         STOP
      END IF

C-----  Implicit integration parameter: THETA
      THETA=PROPS(145)

C-----  Finite deformation ?
C-----  NLGEOM = 0,   small deformation theory
C       otherwise, theory of finite rotation and finite strain, Users
C     must declare "NLGEOM" in the input file, at the *STEP card
C
      IF (PROPS(146).EQ.0.) THEN
         NLGEOM=0
      ELSE
         NLGEOM=1
      END IF

C-----  Iteration?
C-----  ITRATN = 0, no iteration
C       otherwise, iteration (solving increments of stresses and
C     solution dependent state variables)
C
      IF (PROPS(153).EQ.0.) THEN
         ITRATN=0
      ELSE
         ITRATN=1
      END IF

      ITRMAX=NINT(PROPS(154))
      GAMERR=PROPS(155)

      NITRTN=-1

      DO I=1,ndir+nshr
         DSOLD(I)=0.
      END DO

      DO J=1,ND
         DGAMOD(J)=0.
         DTAUOD(J)=0.
         DGSPOD(J)=0.
         DO I=1,3
            DSPNRO(I,J)=0.
            DSPDRO(I,J)=0.
         END DO
      END DO

C-----  Increment of spin associated with the material element: DSPIN
C     (only needed for finite rotation)
C
      IF (NLGEOM.NE.0) THEN

         DSPIN(1)=0.0
         DSPIN(2)=0.0
         DSPIN(3)=0.0

      END IF

C-----  Increment of dilatational strain: DEV

      DEV=0.D0
      DO I=1,ndir
         DEV=DEV+DSTRAN(I)
      END DO

C-----  Iteration starts (only when iteration method is used)

1000  CONTINUE

C-----  Parameter NITRTN: number of iterations
C       NITRTN = 0 --- no-iteration solution
C
      NITRTN=NITRTN+1

C-----  Check whether the current stress state is the initial state

      IF (STATEV(1).EQ.0.) THEN

C-----  Initial state
C
C-----  Generating the following parameters and variables at initial
C     state:
C          Total number of slip systems in all the sets NSLPTL
C          Number of slip systems in each set NSLIP
C          Unit vectors in initial slip directions SLPDIR
C          Unit normals to initial slip planes SLPNOR
C
         NSLPTL=0
         DO I=1,NSET
            ISPNOR(1)=NINT(PROPS(25+8*I))
            ISPNOR(2)=NINT(PROPS(26+8*I))
            ISPNOR(3)=NINT(PROPS(27+8*I))

            ISPDIR(1)=NINT(PROPS(28+8*I))
            ISPDIR(2)=NINT(PROPS(29+8*I))
            ISPDIR(3)=NINT(PROPS(30+8*I))

            CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),
     2                    SLPNOR(1,NSLPTL+1), ROTATE)

            NSLPTL=NSLPTL+NSLIP(I)
         END DO

         IF (ND.LT.NSLPTL) THEN
            WRITE (6,*)
     2 '***ERROR - parameter ND chosen by the present user is less than
     3             the total number of slip systems NSLPTL'
            STOP
         END IF

C-----  Slip deformation tensor: SLPDEF (Schmid factors)
         DO J=1,NSLPTL
            SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)
            SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)
            SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)
            SLPDEF(4,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)
     .                        +SLPDIR(2,J)*SLPNOR(1,J))
            SLPDEF(6,J)=0.5*(SLPDIR(1,J)*SLPNOR(3,J)
     .                        +SLPDIR(3,J)*SLPNOR(1,J))
            SLPDEF(5,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)
     .                        +SLPDIR(3,J)*SLPNOR(2,J))
         END DO

C-----  Initial value of state variables: unit normal to a slip plane
C     and unit vector in a slip direction
C
C     ---change ---nstatev from 125 to 135 
         STATEV(nstatev-10)=FLOAT(NSLPTL)
         DO I=1,NSET
C     ---change ---nstatev from 125 to 135  #########
            STATEV(nstatev-10-4+I)=FLOAT(NSLIP(I))
         END DO

         IDNOR=3*NSLPTL
         IDDIR=6*NSLPTL
         DO J=1,NSLPTL
            DO I=1,3
               IDNOR=IDNOR+1
               STATEV(IDNOR)=SLPNOR(I,J)

               IDDIR=IDDIR+1
               STATEV(IDDIR)=SLPDIR(I,J)
            END DO
         END DO

C-----  Initial value of the current strength for all slip systems
C
         CALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))

C-----  Initial value of shear strain in slip systems
CFIX--  Initial value of cumulative shear strain in each slip systems

         DO I=1,NSLPTL
            STATEV(NSLPTL+I)=0.
CFIXA
            STATEV(9*NSLPTL+I)=0.
CFIXB
         END DO

CFIXA
         STATEV(10*NSLPTL+1)=0.
CFIXB
   
C     ---change ---nstatev from 125 to 135 
C--------Initial value for state variable 126 to 135
    DO I=1,10
    STATEV(125+I)=0.0
    END DO

C-----  Initial value of the resolved shear stress in slip systems
         DO I=1,NSLPTL
            TERM1=0.

            DO J=1,ndir+nshr
               IF (J.LE.ndir) THEN
                  TERM1=TERM1+SLPDEF(J,I)*STRESS(J)
               ELSE
                  TERM1=TERM1+SLPDEF(J-ndir+3,I)*STRESS(J)
               END IF
            END DO

            STATEV(2*NSLPTL+I)=TERM1
         END DO

      ELSE

C-----  Current stress state
C
C-----  Copying from the array of state variables STATVE the following
C          parameters and variables at current stress state:
C          Total number of slip systems in all the sets NSLPTL
C          Number of slip systems in each set NSLIP
C          Current slip directions SLPDIR
C          Normals to current slip planes SLPNOR
C
C     ---change ---nstatev from 125 to 135  #####
         NSLPTL=NINT(STATEV(nstatev-10))
         DO I=1,NSET
            NSLIP(I)=NINT(STATEV(nstatev-14+I))
         END DO

         IDNOR=3*NSLPTL
         IDDIR=6*NSLPTL
         DO J=1,NSLPTL
            DO I=1,3
               IDNOR=IDNOR+1
               SLPNOR(I,J)=STATEV(IDNOR)

               IDDIR=IDDIR+1
               SLPDIR(I,J)=STATEV(IDDIR)
            END DO
         END DO

C-----  Slip deformation tensor: SLPDEF (Schmid factors)
         DO J=1,NSLPTL
            SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)
            SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)
            SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)
            SLPDEF(4,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)
     .                    +SLPDIR(2,J)*SLPNOR(1,J))
            SLPDEF(6,J)=0.5*(SLPDIR(1,J)*SLPNOR(3,J)
     .                    +SLPDIR(3,J)*SLPNOR(1,J))
            SLPDEF(5,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)
     .                    +SLPDIR(3,J)*SLPNOR(2,J))
         END DO

      END IF

C-----  Slip spin tensor: SLPSPN (only needed for finite rotation)
      IF (NLGEOM.NE.0) THEN
         DO J=1,NSLPTL
            SLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-
     2                       SLPDIR(2,J)*SLPNOR(1,J))
            SLPSPN(3,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-
     2                       SLPDIR(1,J)*SLPNOR(3,J))
            SLPSPN(2,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-
     2                       SLPDIR(3,J)*SLPNOR(2,J))
         END DO
      END IF

C-----  Double dot product of elastic moduli tensor with the slip
C     deformation tensor (Schmid factors) plus, only for finite
C     rotation, the dot product of slip spin tensor with the stress:
C     DDEMSD
C
      DO J=1,NSLPTL
         DO I=1,6
            DDEMSD(I,J)=0.
            DO K=1,6
               DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)
            END DO
         END DO
      END DO

      IF (NLGEOM.NE.0) THEN

        DO J=1,NSLPTL

        DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*stress(1)
        DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(3,J)*stress(1)

        IF (ndir.GT.1) THEN
        DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*stress(2)
        DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(2,J)*stress(2)
        END IF

        IF (ndir.GT.2) THEN
        DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*stress(3)
        DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(3,J)*stress(3)
        END IF

        IF (nshr.GE.1) THEN
        DDEMSD(1,J)=DDEMSD(1,J)+2.*SLPSPN(1,J)*stress(ndir+1)
        DDEMSD(2,J)=DDEMSD(2,J)-2.*SLPSPN(1,J)*stress(ndir+1)
        DDEMSD(5,J)=DDEMSD(5,J)+   SLPSPN(3,J)*stress(ndir+1)
        DDEMSD(6,J)=DDEMSD(6,J)-   SLPSPN(2,J)*stress(ndir+1)
        END IF

        IF (nshr.GE.2) THEN
        DDEMSD(2,J)=DDEMSD(2,J)+2.*SLPSPN(2,J)*stress(ndir+2)
        DDEMSD(3,J)=DDEMSD(3,J)-2.*SLPSPN(2,J)*stress(ndir+2)
        DDEMSD(4,J)=DDEMSD(4,J)-   SLPSPN(3,J)*stress(ndir+2)
        DDEMSD(6,J)=DDEMSD(6,J)+   SLPSPN(1,J)*stress(ndir+2)
        END IF

        IF (nshr.EQ.3) THEN
        DDEMSD(1,J)=DDEMSD(1,J)-2.*SLPSPN(3,J)*stress(ndir+3)
        DDEMSD(3,J)=DDEMSD(3,J)+2.*SLPSPN(3,J)*stress(ndir+3)
        DDEMSD(4,J)=DDEMSD(4,J)+   SLPSPN(2,J)*stress(ndir+3)
        DDEMSD(5,J)=DDEMSD(5,J)-   SLPSPN(1,J)*stress(ndir+3)
        END IF
       
      END DO

      END IF

C-----  Shear strain-rate in a slip system at the start of increment:
C     FSLIP, and its derivative: DFDXSP
C
      ID=1
      DO I=1,NSET
         IF (I.GT.1) ID=ID+NSLIP(I-1)
         CALL STRAINRATE (STATEV(NSLPTL+ID), STATEV(2*NSLPTL+ID),
     2                    STATEV(ID), NSLIP(I), FSLIP(ID), DFDXSP(ID),
     3                    PROPS(65+8*I))
      END DO

C-----  Self- and latent-hardening laws
CFIXA 
       CALL LATENTHARDEN (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1),
     2                   STATEV(1), STATEV(9*NSLPTL+1),
     3                   STATEV(10*NSLPTL+1), NSLIP, NSLPTL,
     4                   NSET, H(1,1), PROPS(97), ND)
CFIXB

C-----  LU decomposition to solve the increment of shear strain in a
C     slip system
C
      TERM1=THETA*dt
      DO I=1,NSLPTL
         TAUSLP=STATEV(2*NSLPTL+I)
         GSLIP=STATEV(I)
         X=TAUSLP/GSLIP
         TERM2=TERM1*DFDXSP(I)/GSLIP
         TERM3=TERM1*X*DFDXSP(I)/GSLIP

         DO J=1,NSLPTL
            TERM4=0.
            DO K=1,6
               TERM4=TERM4+DDEMSD(K,I)*SLPDEF(K,J)
            END DO

            WORKST(I,J)=TERM2*TERM4+H(I,J)*TERM3*SIGN(1.0,FSLIP(J))

            IF (NITRTN.GT.0) WORKST(I,J)=WORKST(I,J)+TERM3*DHDGDG(I,J)

         END DO

         WORKST(I,I)=WORKST(I,I)+1.
      END DO

      CALL LUDCMP (WORKST, NSLPTL, ND, INDX, DDCMP)

C-----  Increment of shear strain in a slip system: DGAMMA

      TERM1=THETA*dt
      DO I=1,NSLPTL

         IF (NITRTN.EQ.0) THEN
            TAUSLP=STATEV(2*NSLPTL+I)
            GSLIP=STATEV(I)
            X=TAUSLP/GSLIP
            TERM2=TERM1*DFDXSP(I)/GSLIP

            DGAMMA(I)=0.
            DO J=1,ndir
               DGAMMA(I)=DGAMMA(I)+DDEMSD(J,I)*DSTRAN(J)
            END DO

            IF (NSHR.GT.0) THEN
               DO J=1,NSHR
                  DGAMMA(I)=DGAMMA(I)+DDEMSD(J+3,I)*DSTRAN(J+ndir)
               END DO
            END IF

            DGAMMA(I)=DGAMMA(I)*TERM2+FSLIP(I)*dt

         ELSE
            DGAMMA(I)=TERM1*(FSLIP(I)-FSLIP1(I))+FSLIP1(I)*dt
     2                -DGAMOD(I)

         END IF

      END DO

      CALL LUBKSB (WORKST, NSLPTL, ND, INDX, DGAMMA)

      DO I=1,NSLPTL
         DGAMMA(I)=DGAMMA(I)+DGAMOD(I)
      END DO

C-----  Update the shear strain in a slip system: STATEV(NSLPTL+1) -
C     STATEV(2*NSLPTL)
C
      DO I=1,NSLPTL
         STATEV(NSLPTL+I)=STATEV(NSLPTL+I)+DGAMMA(I)-DGAMOD(I)
      END DO

C-----  Increment of current strength in a slip system: DGSLIP
      DO I=1,NSLPTL
         DGSLIP(I)=0.
         DO J=1,NSLPTL
            DGSLIP(I)=DGSLIP(I)+H(I,J)*ABS(DGAMMA(J))
         END DO
      END DO

C-----  Update the current strength in a slip system: STATEV(1) -
C     STATEV(NSLPTL)
C
      DO I=1,NSLPTL
         STATEV(I)=STATEV(I)+DGSLIP(I)-DGSPOD(I)
      END DO

C-----  Increment of strain associated with lattice stretching: DELATS
      DO J=1,6
         DELATS(J)=0.
      END DO

      DO J=1,3
         IF (J.LE.ndir) DELATS(J)=DSTRAN(J)
         DO I=1,NSLPTL
            DELATS(J)=DELATS(J)-SLPDEF(J,I)*DGAMMA(I)
         END DO
      END DO

      DO J=1,3
         IF (J.LE.NSHR) DELATS(J+3)=DSTRAN(J+ndir)
         DO I=1,NSLPTL
            DELATS(J+3)=DELATS(J+3)-SLPDEF(J+3,I)*DGAMMA(I)
         END DO
      END DO

C-----  Increment of deformation gradient associated with lattice
C     stretching in the current state, i.e. the velocity gradient
C     (associated with lattice stretching) times the increment of time:
C     DVGRAD (only needed for finite rotation)
C
      IF (NLGEOM.NE.0) THEN
         DO J=1,3
            DO I=1,3
               IF (I.EQ.J) THEN
                  DVGRAD(I,J)=DELATS(I)
               ELSE
                  DVGRAD(I,J)=DELATS(I+J+1)
               END IF
            END DO
         END DO

         DO J=1,3
            DO I=1,J
               IF (J.GT.I) THEN
                  IJ2=I+J-2
                  IF (MOD(IJ2,2).EQ.1) THEN
                     TERM1=1.
                  ELSE
                     TERM1=-1.
                  END IF

                  DVGRAD(I,J)=DVGRAD(I,J)+TERM1*DSPIN(IJ2)
                  DVGRAD(J,I)=DVGRAD(J,I)-TERM1*DSPIN(IJ2)

                  DO K=1,NSLPTL
                     DVGRAD(I,J)=DVGRAD(I,J)-TERM1*DGAMMA(K)*
     2                                       SLPSPN(IJ2,K)
                     DVGRAD(J,I)=DVGRAD(J,I)+TERM1*DGAMMA(K)*
     2                                       SLPSPN(IJ2,K)
                  END DO
               END IF

            END DO
         END DO

      END IF

C-----  Increment of resolved shear stress in a slip system: DTAUSP
      DO I=1,NSLPTL
         DTAUSP(I)=0.
         DO J=1,6
            DTAUSP(I)=DTAUSP(I)+DDEMSD(J,I)*DELATS(J)
         END DO
      END DO

C-----  Update the resolved shear stress in a slip system:
C     STATEV(2*NSLPTL+1) - STATEV(3*NSLPTL)
C
      DO I=1,NSLPTL
         STATEV(2*NSLPTL+I)=STATEV(2*NSLPTL+I)+DTAUSP(I)-DTAUOD(I)
      END DO

C-----  Increment of stress: DSTRES
      IF (NLGEOM.EQ.0) THEN
         DO I=1,ndir+nshr
            DSTRES(I)=0.
         END DO
      ELSE
         DO I=1,ndir+nshr
            DSTRES(I)=-STRESS(I)*DEV
         END DO
      END IF

      DO I=1,ndir
         DO J=1,ndir
            DSTRES(I)=DSTRES(I)+D(I,J)*DSTRAN(J)
         END DO

         IF (NSHR.GT.0) THEN
            DO J=1,NSHR
               DSTRES(I)=DSTRES(I)+D(I,J+3)*DSTRAN(J+ndir)
            END DO
         END IF

         DO J=1,NSLPTL
            DSTRES(I)=DSTRES(I)-DDEMSD(I,J)*DGAMMA(J)
         END DO
      END DO

      IF (NSHR.GT.0) THEN
         DO I=1,NSHR

            DO J=1,ndir
               DSTRES(I+ndir)=DSTRES(I+ndir)+D(I+3,J)*DSTRAN(J)
            END DO

            DO J=1,NSHR
               DSTRES(I+ndir)=DSTRES(I+ndir)+D(I+3,J+3)*DSTRAN(J+ndir)
            END DO

            DO J=1,NSLPTL
               DSTRES(I+ndir)=DSTRES(I+ndir)-DDEMSD(I+3,J)*DGAMMA(J)
            END DO

         END DO
      END IF
c
C-----  Update the stress: STRESS
c
      IF (steptime.eq.0.) then
    twomu=props(1)/(1.0+props(2))
    alamda=twomu*props(2)/(1.0-2.0*props(2))
    trace=dstran(1)+dstran(2)+dstran(3)
         DO I=1,ndir
            stress(I)=twomu*dstran(I)+alamda*trace
         END DO
        
         DO J=ndir+1, ndir+nshr
            stress(I)=twomu*dstran(I)
         END DO
            
      ELSE

        DO I=1,ndir+nshr
          stress(I)=stress(I)+DSTRES(I)-DSOLD(I)
        END DO
      END IF

c
C-----  Increment of normal to a slip plane and a slip direction (only
C     needed for finite rotation)
C
      IF (NLGEOM.NE.0) THEN
         DO J=1,NSLPTL
            DO I=1,3
               DSPNOR(I,J)=0.
               DSPDIR(I,J)=0.

               DO K=1,3
                  DSPNOR(I,J)=DSPNOR(I,J)-SLPNOR(K,J)*DVGRAD(K,I)
                  DSPDIR(I,J)=DSPDIR(I,J)+SLPDIR(K,J)*DVGRAD(I,K)
               END DO

            END DO
         END DO

C-----  Update the normal to a slip plane and a slip direction (only
C     needed for finite rotation)
C
         IDNOR=3*NSLPTL
         IDDIR=6*NSLPTL
         DO J=1,NSLPTL
            DO I=1,3
               IDNOR=IDNOR+1
               STATEV(IDNOR)=STATEV(IDNOR)+DSPNOR(I,J)-DSPNRO(I,J)

               IDDIR=IDDIR+1
               STATEV(IDDIR)=STATEV(IDDIR)+DSPDIR(I,J)-DSPDRO(I,J)
            END DO
         END DO

      END IF
c
C-----  Iteration ?
c
      IF (ITRATN.NE.0) THEN

C-----  Save solutions (without iteration):
C            Shear strain-rate in a slip system FSLIP1
C            Current strength in a slip system GSLP1
C            Shear strain in a slip system GAMMA1
C            Resolved shear stress in a slip system TAUSP1
C            Normal to a slip plane SPNOR1
C            Slip direction SPDIR1
C            Stress STRES1
C            Jacobian matrix DDSDE1
C
         IF (NITRTN.EQ.0) THEN

            IDNOR=3*NSLPTL
            IDDIR=6*NSLPTL
            DO J=1,NSLPTL
               FSLIP1(J)=FSLIP(J)
               GSLP1(J)=STATEV(J)
               GAMMA1(J)=STATEV(NSLPTL+J)
               TAUSP1(J)=STATEV(2*NSLPTL+J)
               DO I=1,3
                  IDNOR=IDNOR+1
                  SPNOR1(I,J)=STATEV(IDNOR)

                  IDDIR=IDDIR+1
                  SPDIR1(I,J)=STATEV(IDDIR)
               END DO
            END DO

            DO J=1,ndir+nshr
               STRES1(J)=STRESS(J)
cc               DO I=1,NTENS
cc                  DDSDE1(I,J)=DDSDDE(I,J)
cc               END DO
            END DO

         END IF

C-----  Increments of stress DSOLD, and solution dependent state
C     variables DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO (for the next
C     iteration)
C
         DO I=1,ndir+nshr
            DSOLD(I)=DSTRES(I)
         END DO

         DO J=1,NSLPTL
            DGAMOD(J)=DGAMMA(J)
            DTAUOD(J)=DTAUSP(J)
            DGSPOD(J)=DGSLIP(J)
            DO I=1,3
               DSPNRO(I,J)=DSPNOR(I,J)
               DSPDRO(I,J)=DSPDIR(I,J)
            END DO
         END DO

C-----  Check if the iteration solution converges
         IDBACK=0
         ID=0
         DO I=1,NSET
            DO J=1,NSLIP(I)
               ID=ID+1
               X=STATEV(2*NSLPTL+ID)/STATEV(ID)
               RESIDU=THETA*dt*F(X,PROPS(65+8*I))+dt*(1.0-THETA)*
     2                FSLIP1(ID)-DGAMMA(ID)
               IF (ABS(RESIDU).GT.GAMERR) IDBACK=1
            END DO
         END DO

         IF (IDBACK.NE.0.AND.NITRTN.LT.ITRMAX) THEN
C-----  Iteration: arrays for iteration
CFIXA
            CALL ITERATION (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1),
     2                      STATEV(1), STATEV(9*NSLPTL+1),
     3                      STATEV(10*NSLPTL+1), NSLPTL,
     4                      NSET, NSLIP, ND, PROPS(97), DGAMOD,
     5                      DHDGDG)
CFIXB

            GO TO 1000

         ELSE IF (NITRTN.GE.ITRMAX) THEN
C-----  Solution not converge within maximum number of iteration (the
C     solution without iteration will be used)
C
            DO J=1,ndir+nshr
               STRESS(J)=STRES1(J)
            END DO

            IDNOR=3*NSLPTL
            IDDIR=6*NSLPTL
            DO J=1,NSLPTL
               STATEV(J)=GSLP1(J)
               STATEV(NSLPTL+J)=GAMMA1(J)
               STATEV(2*NSLPTL+J)=TAUSP1(J)

               DO I=1,3
                  IDNOR=IDNOR+1
                  STATEV(IDNOR)=SPNOR1(I,J)

                  IDDIR=IDDIR+1
                  STATEV(IDDIR)=SPDIR1(I,J)
               END DO
            END DO

         END IF

      END IF

c
CFIX--  Total cumulative shear strains on each slip system (sum of the
CFIX    absolute values of shear strains in each individual slip system)
C
      DO I=1,NSLPTL
CFIXA
         STATEV(10*NSLPTL+1)=STATEV(10*NSLPTL+1)+ABS(DGAMMA(I))
        
         STATEV(9*NSLPTL+I)=STATEV(9*NSLPTL+I)+ABS(DGAMMA(I))
CFIXB
      END DO

        STATEV(NSTATEV-5)=STATEV(10*NSLPTL+1)
c-------store stressNew, stateNew --------------

      DO NK8=1,ndir+nshr
        STRESSNEW(NK9,NK8)=STRESS(NK8)
      END DO

      DO NK8=1,nstatev
        STATENEW(NK9,NK8)=STATEV(NK8)
      END DO
c
c    update the specific internal energy
c
    stressPower=0.0
    do i=1, ndir
        stressPower=stressPower+0.5*(strainInc(nk9,i)*
     .                (stressOld(nk9,i)+stressNew(nk9,i)))
      end do
    do i=1, nshr
    j=i+nshr
        stressPower=stressPower+(stressOld(nk9,j)+stressNew(nk9,j))*
     .                strainInc(nk9,j)
    end do
c
    enerInternNew(nk9)=enerInternOld(nk9)+stressPower/density(nk9)
c
c    update the dissipated inelastic specific energy
c
    smean=(stressNew(nk9,1)+stressNew(nk9,2)+stressNew(nk9,3))/3.0
    equivStress=0.0
    do i=1, ndir
        equivStress=equivStress+(stressNew(nk9,i)-smean)**2
    end do
    do i=1, nshr
        equivStress=equivStress+2.0*stressNew(nk9,i+ndir)**2
    end do
    equivStress=sqrt(1.5*equivStress)
c
c schmit factor
c
      IDNOR=3*NSLPTL
      IDDIR=6*NSLPTL
      DO J=1,NSLPTL
    DGAMMA(j)=stateNew(nk9,j+nslptl)-stateOld(nk9,j+nslptl)
         DO I=1,3
            IDNOR=IDNOR+1
            SLPNOR(I,J)=stateNew(nk9,IDNOR)
            IDDIR=IDDIR+1
            SLPDIR(I,J)=stateNew(nk9,IDDIR)
         END DO
      END DO

C-----  Slip deformation tensor: SLPDEF (Schmid factors)

      DO J=1,NSLPTL
         SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)
         SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)
         SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)
         SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)
         SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)
         SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J)
      END DO
c
      DO J=1,6
         DELATS(J)=0.
      END DO
c
      DO J=1,3
         DO I=1,NSLPTL
            DELATS(J)=DELATS(J)+SLPDEF(J,I)*DGAMMA(I)
         END DO
      END DO
c
      DO J=1,3
         DO I=1,NSLPTL
            DELATS(J+3)=DELATS(J+3)+SLPDEF(J+3,I)*DGAMMA(I)
         END DO
      END DO

    deqps=sqrt(2*(delats(1)**2+delats(2)**2+delats(3)**2+
     .              2*(delats(4)**2+delats(5)**2+delats(6)**2))/3.0)
    plasticWorkInc=equivStress*deqps
    enerInelasNew(nk9)=enerInelasOld(nk9)+plasticWorkInc/density(nk9)
c
100    CONTINUE

      RETURN
      END

C----------------------------------------------------------------------

      SUBROUTINE ROTATION (PROP, ROTATE)

C-----  This subroutine calculates the rotation matrix, i.e. the
C     direction cosines of cubic crystal [100], [010] and [001]
C     directions in global system

C-----  The rotation matrix is stored in the array ROTATE.

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
CC      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION PROP(16), ROTATE(3,3), TERM1(3,3), TERM2(3,3), INDX(3)

C-----  Subroutines:
C
C       CROSS  -- cross product of two vectors
C
C       LUDCMP -- LU decomposition
C
C       LUBKSB -- linear equation solver based on LU decomposition
C                 method (must call LUDCMP first)

C-----  PROP -- constants characterizing the crystal orientation
C               (INPUT)
C
C            PROP(1) - PROP(3) -- direction of the first vector in
C                                 local cubic crystal system
C            PROP(4) - PROP(6) -- direction of the first vector in
C                                 global system
C
C            PROP(9) - PROP(11)-- direction of the second vector in
C                                 local cubic crystal system
C            PROP(12)- PROP(14)-- direction of the second vector in
C                                 global system
C
C-----  ROTATE -- rotation matrix (OUTPUT):
C
C            ROTATE(i,1) -- direction cosines of direction [1 0 0] in
C                           local cubic crystal system
C            ROTATE(i,2) -- direction cosines of direction [0 1 0] in
C                           local cubic crystal system
C            ROTATE(i,3) -- direction cosines of direction [0 0 1] in
C                           local cubic crystal system

C-----  local matrix: TERM1
      CALL CROSS (PROP(1), PROP(9), TERM1, ANGLE1)

C-----  LU decomposition of TERM1
      CALL LUDCMP (TERM1, 3, 3, INDX, DCMP)

C-----  inverse matrix of TERM1: TERM2
      DO J=1,3
         DO I=1,3
            IF (I.EQ.J) THEN
               TERM2(I,J)=1.
            ELSE
               TERM2(I,J)=0.
            END IF
         END DO
      END DO

      DO J=1,3
         CALL LUBKSB (TERM1, 3, 3, INDX, TERM2(1,J))
      END DO

C-----  global matrix: TERM1
      CALL CROSS (PROP(4), PROP(12), TERM1, ANGLE2)

C-----  Check: the angle between first and second vector in local and
C     global systems must be the same.  The relative difference must be
C     less than 0.1%.
C
      IF (ABS(ANGLE1/ANGLE2-1.).GT.0.001) THEN
         WRITE (6,*)
     2      '***ERROR - angles between two vectors are not the same'
         STOP
      END IF

C-----  rotation matrix: ROTATE
      DO J=1,3
         DO I=1,3
            ROTATE(I,J)=0.
            DO K=1,3
               ROTATE(I,J)=ROTATE(I,J)+TERM1(I,K)*TERM2(K,J)
            END DO
         END DO
      END DO

      RETURN
      END

C-----------------------------------

           SUBROUTINE CROSS (A, B, C, ANGLE)

C-----  (1) normalize vectors A and B to unit vectors
C       (2) store A, B and A*B (cross product) in C

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
           DIMENSION A(3), B(3), C(3,3)

           SUM1=SQRT(A(1)**2+A(2)**2+A(3)**2)
           SUM2=SQRT(B(1)**2+B(2)**2+B(3)**2)

           IF (SUM1.EQ.0.) THEN
              WRITE (6,*) '***ERROR - first vector is zero'
              STOP
           ELSE
              DO I=1,3
                 C(I,1)=A(I)/SUM1
              END DO
           END IF

           IF (SUM2.EQ.0.) THEN
              WRITE (6,*) '***ERROR - second vector is zero'
              STOP
           ELSE
              DO I=1,3
                 C(I,2)=B(I)/SUM2
              END DO
           END IF

           ANGLE=0.
           DO I=1,3
              ANGLE=ANGLE+C(I,1)*C(I,2)
           END DO
           ANGLE=ACOS(ANGLE)

           C(1,3)=C(2,1)*C(3,2)-C(3,1)*C(2,2)
           C(2,3)=C(3,1)*C(1,2)-C(1,1)*C(3,2)
           C(3,3)=C(1,1)*C(2,2)-C(2,1)*C(1,2)
           SUM3=SQRT(C(1,3)**2+C(2,3)**2+C(3,3)**2)
           IF (SUM3.LT.1.E-8) THEN
              WRITE (6,*)
     2           '***ERROR - first and second vectors are parallel'
               STOP
            END IF

           RETURN
           END

C----------------------------------------------------------------------

      SUBROUTINE SLIPSYS (ISPDIR, ISPNOR, NSLIP, SLPDIR, SLPNOR,
     2                    ROTATE)

C-----  This subroutine generates all slip systems in the same set for
C     a CUBIC crystal.  For other crystals (e.g., HCP, Tetragonal,
C     Orthotropic, ...), it has to be modified to include the effect of
C     crystal aspect ratio.

C-----  Denote s as a slip direction and m as normal to a slip plane. 
C     In a cubic crystal, (s,-m), (-s,m) and (-s,-m) are NOT considered
C     independent of (s,m).

C-----  Subroutines:  LINE1 and LINE

C-----  Variables:
C
C     ISPDIR -- a typical slip direction in this set of slip systems
C               (integer)  (INPUT)
C     ISPNOR -- a typical normal to slip plane in this set of slip
C               systems (integer)  (INPUT)
C     NSLIP  -- number of independent slip systems in this set
C               (OUTPUT)
C     SLPDIR -- unit vectors of all slip directions  (OUTPUT)
C     SLPNOR -- unit normals to all slip planes  (OUTPUT)
C     ROTATE -- rotation matrix (INPUT)
C          ROTATE(i,1) -- direction cosines of [100] in global system
C          ROTATE(i,2) -- direction cosines of [010] in global system
C          ROTATE(i,3) -- direction cosines of [001] in global system
C
C     NSPDIR -- number of all possible slip directions in this set
C     NSPNOR -- number of all possible slip planes in this set
C     IWKDIR -- all possible slip directions (integer)
C     IWKNOR -- all possible slip planes (integer)

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
CC      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ISPDIR(3), ISPNOR(3), SLPDIR(3,50), SLPNOR(3,50),
     *          ROTATE(3,3), IWKDIR(3,24), IWKNOR(3,24), TERM(3)

      NSLIP=0
      NSPDIR=0
      NSPNOR=0

C-----  Generating all possible slip directions in this set
C
C       Denote the slip direction by [lmn].  I1 is the minimum of the
C     absolute value of l, m and n, I3 is the maximum and I2 is the
C     mode, e.g. (1 -3 2), I1=1, I2=2 and I3=3.  I1<=I2<=I3.

      I1=MIN(IABS(ISPDIR(1)),IABS(ISPDIR(2)),IABS(ISPDIR(3)))
      I3=MAX(IABS(ISPDIR(1)),IABS(ISPDIR(2)),IABS(ISPDIR(3)))
      I2=IABS(ISPDIR(1))+IABS(ISPDIR(2))+IABS(ISPDIR(3))-I1-I3

      RMODIR=SQRT(FLOAT(I1*I1+I2*I2+I3*I3))

C     I1=I2=I3=0
      IF (I3.EQ.0) THEN
         WRITE (6,*) '***ERROR - slip direction is [000]'
         STOP

C     I1=I2=0, I3>0   ---   [001] type
      ELSE IF (I2.EQ.0) THEN
         NSPDIR=3
         DO J=1,3
            DO I=1,3
               IWKDIR(I,J)=0
               IF (I.EQ.J) IWKDIR(I,J)=I3
            END DO
         END DO

C     I1=0, I3>=I2>0
      ELSE IF (I1.EQ.0) THEN

C        I1=0, I3=I2>0   ---   [011] type
         IF (I2.EQ.I3) THEN
            NSPDIR=6
            DO J=1,6
               DO I=1,3
                  IWKDIR(I,J)=I2
                  IF (I.EQ.J.OR.J-I.EQ.3) IWKDIR(I,J)=0
                  IWKDIR(1,6)=-I2
                  IWKDIR(2,4)=-I2
                  IWKDIR(3,5)=-I2
               END DO
            END DO

C        I1=0, I3>I2>0   ---   [012] type
         ELSE
            NSPDIR=12
            CALL LINE1 (I2, I3, IWKDIR(1,1), 1)
            CALL LINE1 (I3, I2, IWKDIR(1,3), 1)
            CALL LINE1 (I2, I3, IWKDIR(1,5), 2)
            CALL LINE1 (I3, I2, IWKDIR(1,7), 2)
            CALL LINE1 (I2, I3, IWKDIR(1,9), 3)
            CALL LINE1 (I3, I2, IWKDIR(1,11), 3)

         END IF

C     I1=I2=I3>0   ---   [111] type
      ELSE IF (I1.EQ.I3) THEN
         NSPDIR=4
         CALL LINE (I1, I1, I1, IWKDIR)

C     I3>I2=I1>0   ---   [112] type
      ELSE IF (I1.EQ.I2) THEN
         NSPDIR=12
         CALL LINE (I1, I1, I3, IWKDIR(1,1))
         CALL LINE (I1, I3, I1, IWKDIR(1,5))
         CALL LINE (I3, I1, I1, IWKDIR(1,9))

C     I3=I2>I1>0   ---   [122] type
      ELSE IF (I2.EQ.I3) THEN
         NSPDIR=12
         CALL LINE (I1, I2, I2, IWKDIR(1,1))
         CALL LINE (I2, I1, I2, IWKDIR(1,5))
         CALL LINE (I2, I2, I1, IWKDIR(1,9))

C     I3>I2>I1>0   ---   [123] type
      ELSE
         NSPDIR=24
         CALL LINE (I1, I2, I3, IWKDIR(1,1))
         CALL LINE (I3, I1, I2, IWKDIR(1,5))
         CALL LINE (I2, I3, I1, IWKDIR(1,9))
         CALL LINE (I1, I3, I2, IWKDIR(1,13))
         CALL LINE (I2, I1, I3, IWKDIR(1,17))
         CALL LINE (I3, I2, I1, IWKDIR(1,21))

      END IF

C-----  Generating all possible slip planes in this set
C
C       Denote the normal to slip plane by (pqr).  J1 is the minimum of
C     the absolute value of p, q and r, J3 is the maximum and J2 is the
C     mode, e.g. (1 -2 1), J1=1, J2=1 and J3=2.  J1<=J2<=J3.

      J1=MIN(IABS(ISPNOR(1)),IABS(ISPNOR(2)),IABS(ISPNOR(3)))
      J3=MAX(IABS(ISPNOR(1)),IABS(ISPNOR(2)),IABS(ISPNOR(3)))
      J2=IABS(ISPNOR(1))+IABS(ISPNOR(2))+IABS(ISPNOR(3))-J1-J3

      RMONOR=SQRT(FLOAT(J1*J1+J2*J2+J3*J3))

      IF (J3.EQ.0) THEN
         WRITE (6,*) '***ERROR - slip plane is [000]'
         STOP

C     (001) type
      ELSE IF (J2.EQ.0) THEN
         NSPNOR=3
         DO J=1,3
            DO I=1,3
               IWKNOR(I,J)=0
               IF (I.EQ.J) IWKNOR(I,J)=J3
            END DO
         END DO

      ELSE IF (J1.EQ.0) THEN

C     (011) type
         IF (J2.EQ.J3) THEN
            NSPNOR=6
            DO J=1,6
               DO I=1,3
                  IWKNOR(I,J)=J2
                  IF (I.EQ.J.OR.J-I.EQ.3) IWKNOR(I,J)=0
                  IWKNOR(1,6)=-J2
                  IWKNOR(2,4)=-J2
                  IWKNOR(3,5)=-J2
               END DO
            END DO

C     (012) type
         ELSE
            NSPNOR=12
            CALL LINE1 (J2, J3, IWKNOR(1,1), 1)
            CALL LINE1 (J3, J2, IWKNOR(1,3), 1)
            CALL LINE1 (J2, J3, IWKNOR(1,5), 2)
            CALL LINE1 (J3, J2, IWKNOR(1,7), 2)
            CALL LINE1 (J2, J3, IWKNOR(1,9), 3)
            CALL LINE1 (J3, J2, IWKNOR(1,11), 3)

         END IF

C     (111) type
      ELSE IF (J1.EQ.J3) THEN
         NSPNOR=4
         CALL LINE (J1, J1, J1, IWKNOR)

C     (112) type
      ELSE IF (J1.EQ.J2) THEN
         NSPNOR=12
         CALL LINE (J1, J1, J3, IWKNOR(1,1))
         CALL LINE (J1, J3, J1, IWKNOR(1,5))
         CALL LINE (J3, J1, J1, IWKNOR(1,9))

C     (122) type
      ELSE IF (J2.EQ.J3) THEN
         NSPNOR=12
         CALL LINE (J1, J2, J2, IWKNOR(1,1))
         CALL LINE (J2, J1, J2, IWKNOR(1,5))
         CALL LINE (J2, J2, J1, IWKNOR(1,9))

C     (123) type
      ELSE
         NSPNOR=24
         CALL LINE (J1, J2, J3, IWKNOR(1,1))
         CALL LINE (J3, J1, J2, IWKNOR(1,5))
         CALL LINE (J2, J3, J1, IWKNOR(1,9))
         CALL LINE (J1, J3, J2, IWKNOR(1,13))
         CALL LINE (J2, J1, J3, IWKNOR(1,17))
         CALL LINE (J3, J2, J1, IWKNOR(1,21))

      END IF

C-----  Generating all slip systems in this set
C
C-----  Unit vectors in slip directions: SLPDIR, and unit normals to
C     slip planes: SLPNOR in local cubic crystal system
C
CC      WRITE (6,*) '          '
CC      WRITE (6,*) ' #          Slip plane          Slip direction'

      DO J=1,NSPNOR
         DO I=1,NSPDIR

            IDOT=0
            DO K=1,3
               IDOT=IDOT+IWKDIR(K,I)*IWKNOR(K,J)
            END DO

            IF (IDOT.EQ.0) THEN
               NSLIP=NSLIP+1
               DO K=1,3
                  SLPDIR(K,NSLIP)=IWKDIR(K,I)/RMODIR
                  SLPNOR(K,NSLIP)=IWKNOR(K,J)/RMONOR
               END DO

CC               WRITE (6,10) NSLIP,
CC     2                      (IWKNOR(K,J),K=1,3), (IWKDIR(K,I),K=1,3)

            END IF

         END DO
      END DO
CC10    FORMAT(1X,I2,9X,'(',3(1X,I2),1X,')',10X,'[',3(1X,I2),1X,']')

CC      WRITE (6,*) 'Number of slip systems in this set = ',NSLIP
CC      WRITE (6,*) '          '

      IF (NSLIP.EQ.0) THEN
         WRITE (6,*)
     *      'There is no slip direction normal to the slip planes!'
         STOP

      ELSE

C-----  Unit vectors in slip directions: SLPDIR, and unit normals to
C     slip planes: SLPNOR in global system
C
         DO J=1,NSLIP
            DO I=1,3
               TERM(I)=0.
               DO K=1,3
                  TERM(I)=TERM(I)+ROTATE(I,K)*SLPDIR(K,J)
               END DO
            END DO
            DO I=1,3
               SLPDIR(I,J)=TERM(I)
            END DO

            DO I=1,3
               TERM(I)=0.
               DO K=1,3
                  TERM(I)=TERM(I)+ROTATE(I,K)*SLPNOR(K,J)
               END DO
            END DO
            DO I=1,3
               SLPNOR(I,J)=TERM(I)
            END DO
         END DO

      END IF

      RETURN
      END

C----------------------------------

           SUBROUTINE LINE (I1, I2, I3, IARRAY)

C-----  Generating all possible slip directions <lmn> (or slip planes
C     {lmn}) for a cubic crystal, where l,m,n are not zeros.

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
           DIMENSION IARRAY(3,4)

           DO J=1,4
              IARRAY(1,J)=I1
              IARRAY(2,J)=I2
              IARRAY(3,J)=I3
           END DO

           DO I=1,3
              DO J=1,4
                 IF (J.EQ.I+1) IARRAY(I,J)=-IARRAY(I,J)
              END DO
           END DO

           RETURN
           END

C-----------------------------------

           SUBROUTINE LINE1 (J1, J2, IARRAY, ID)

C-----  Generating all possible slip directions <0mn> (or slip planes
C     {0mn}) for a cubic crystal, where m,n are not zeros and m does
C     not equal n.

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
           DIMENSION IARRAY(3,2)

           IARRAY(ID,1)=0
           IARRAY(ID,2)=0

           ID1=ID+1
           IF (ID1.GT.3) ID1=ID1-3
           IARRAY(ID1,1)=J1
           IARRAY(ID1,2)=J1

           ID2=ID+2
           IF (ID2.GT.3) ID2=ID2-3
           IARRAY(ID2,1)=J2
           IARRAY(ID2,2)=-J2
 
           RETURN
           END

C----------------------------------------------------------------------

      SUBROUTINE GSLPINIT (GSLIP0, NSLIP, NSLPTL, NSET, PROP)

C-----  This subroutine calculates the initial value of current
C     strength for each slip system in a rate-dependent single crystal.
C     Two sets of initial values, proposed by Asaro, Pierce et al, and
C     by Bassani, respectively, are used here.  Both sets assume that
C     the initial values for all slip systems are the same (initially
C     isotropic).

C-----  These initial values are assumed the same for all slip systems
C     in each set, though they could be different from set to set, e.g.
C     <110>{111} and <110>{100}.

C-----  Users who want to use their own initial values may change the
C     function subprogram GSLP0.  The parameters characterizing these
C     initial values are passed into GSLP0 through array PROP.

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
CC      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL GSLP0
      DIMENSION GSLIP0(NSLPTL), NSLIP(NSET), PROP(16,NSET)

C-----  Function subprograms:
C
C       GSLP0 -- User-supplied function subprogram given the initial
C                value of current strength at initial state

C-----  Variables:
C
C     GSLIP0 -- initial value of current strength (OUTPUT)
C
C     NSLIP  -- number of slip systems in each set (INPUT)
C     NSLPTL -- total number of slip systems in all the sets (INPUT)
C     NSET   -- number of sets of slip systems (INPUT)
C
C     PROP   -- material constants characterizing the initial value of
C               current strength (INPUT)
C
C               For Asaro, Pierce et al's law
C               PROP(1,i) -- initial hardening modulus H0 in the ith
C                            set of slip systems
C               PROP(2,i) -- saturation stress TAUs in the ith set of 
C                            slip systems
C               PROP(3,i) -- initial critical resolved shear stress
C                            TAU0 in the ith set of slip systems
C
C               For Bassani's law
C               PROP(1,i) -- initial hardening modulus H0 in the ith
C                            set of slip systems
C               PROP(2,i) -- stage I stress TAUI in the ith set of 
C                            slip systems (or the breakthrough stress
C                            where large plastic flow initiates)
C               PROP(3,i) -- initial critical resolved shear stress
C                            TAU0 in the ith set of slip systems
C

      ID=0
      DO I=1,NSET
         ISET=I
         DO J=1,NSLIP(I)
            ID=ID+1
            GSLIP0(ID)=GSLP0(NSLPTL,NSET,NSLIP,PROP(1,I),ID,ISET)
         END DO
      END DO

      RETURN
      END

C----------------------------------

C-----  Use single precision on cray
C
           FUNCTION GSLP0(NSLPTL,NSET,NSLIP,PROP,ISLIP,ISET)

C-----     User-supplied function subprogram given the initial value of
C        current strength at initial state

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
           DIMENSION NSLIP(NSET), PROP(16)

           GSLP0=PROP(3)

           RETURN
           END

C----------------------------------------------------------------------

      SUBROUTINE STRAINRATE (GAMMA, TAUSLP, GSLIP, NSLIP, FSLIP,
     2                       DFDXSP, PROP)

C-----  This subroutine calculates the shear strain-rate in each slip
C     system for a rate-dependent single crystal.  The POWER LAW
C     relation between shear strain-rate and resolved shear stress
C     proposed by Hutchinson, Pan and Rice, is used here.

C-----  The power law exponents are assumed the same for all slip
C     systems in each set, though they could be different from set to
C     set, e.g. <110>{111} and <110>{100}.  The strain-rate coefficient
C     in front of the power law form are also assumed the same for all
C     slip systems in each set.

C-----  Users who want to use their own constitutive relation may
C     change the function subprograms F and its derivative DFDX,
C     where F is the strain hardening law, dGAMMA/dt = F(X),
C     X=TAUSLP/GSLIP.  The parameters characterizing F are passed into
C     F and DFDX through array PROP.

C-----  Function subprograms:
C
C       F    -- User-supplied function subprogram which gives shear
C               strain-rate for each slip system based on current
C               values of resolved shear stress and current strength
C
C       DFDX -- User-supplied function subprogram dF/dX, where x is the
C               ratio of resolved shear stress over current strength

C-----  Variables:
C
C     GAMMA  -- shear strain in each slip system at the start of time
C               step  (INPUT)
C     TAUSLP -- resolved shear stress in each slip system (INPUT)
C     GSLIP  -- current strength (INPUT)
C     NSLIP  -- number of slip systems in this set (INPUT)
C
C     FSLIP  -- current value of F for each slip system (OUTPUT)
C     DFDXSP -- current value of DFDX for each slip system (OUTPUT)
C
C     PROP   -- material constants characterizing the strain hardening
C               law (INPUT)
C
C               For the current power law strain hardening law
C               PROP(1) -- power law hardening exponent
C               PROP(1) = infinity corresponds to a rate-independent
C               material
C               PROP(2) -- coefficient in front of power law hardening

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
CC      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL F, DFDX
      DIMENSION GAMMA(NSLIP), TAUSLP(NSLIP), GSLIP(NSLIP),
     2          FSLIP(NSLIP), DFDXSP(NSLIP), PROP(8)

      DO I=1,NSLIP
         X=TAUSLP(I)/GSLIP(I)
         FSLIP(I)=F(X,PROP)
         DFDXSP(I)=DFDX(X,PROP)
      END DO

      RETURN
      END

C-----------------------------------

C-----  Use single precision on cray
C
           FUNCTION F(X,PROP)

C-----     User-supplied function subprogram which gives shear
C        strain-rate for each slip system based on current values of
C        resolved shear stress and current strength
C
C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
           DIMENSION PROP(8)

           F=PROP(2)*(ABS(X))**PROP(1)*SIGN(1.0,X)

           RETURN
           END

C-----------------------------------

C-----  Use single precision on cray
C
           FUNCTION DFDX(X,PROP)

C-----     User-supplied function subprogram dF/dX, where x is the
C        ratio of resolved shear stress over current strength

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
           DIMENSION PROP(8)

           DFDX=PROP(1)*PROP(2)*(ABS(X))**(PROP(1)-1.)

           RETURN
           END

C----------------------------------------------------------------------

CFIXA
      SUBROUTINE LATENTHARDEN (GAMMA, TAUSLP, GSLIP, GMSLTL, GAMTOL,
     2                         NSLIP, NSLPTL, NSET, H, PROP, ND)
CFIXB

C-----  This subroutine calculates the current self- and latent-
C     hardening moduli for all slip systems in a rate-dependent single
C     crystal.  Two kinds of hardening law are used here.  The first
C     law, proposed by Asaro, and Pierce et al, assumes a HYPER SECANT
C     relation between self- and latent-hardening moduli and overall
C     shear strain.  The Bauschinger effect has been neglected.  The
C     second is Bassani's hardening law, which gives an explicit
C     expression of slip interactions between slip systems.  The
C     classical three stage hardening for FCC single crystal could be
C     simulated.

C-----  The hardening coefficients are assumed the same for all slip
C     systems in each set, though they could be different from set to
C     set, e.g. <110>{111} and <110>{100}.

C-----  Users who want to use their own self- and latent-hardening law
C     may change the function subprograms HSELF (self hardening) and
C     HLATNT (latent hardening).  The parameters characterizing these
C     hardening laws are passed into HSELF and HLATNT through array
C     PROP.

C-----  Function subprograms:
C
C       HSELF  -- User-supplied self-hardening function in a slip
C                 system
C
C       HLATNT -- User-supplied latent-hardening function

C-----  Variables:
C
C     GAMMA  -- shear strain in all slip systems at the start of time
C               step  (INPUT)
C     TAUSLP -- resolved shear stress in all slip systems (INPUT)
C     GSLIP  -- current strength (INPUT)
CFIX  GMSLTL -- total cumulative shear strains on each individual slip system
CFIX            (INPUT)
C     GAMTOL -- total cumulative shear strains over all slip systems
C               (INPUT)
C     NSLIP  -- number of slip systems in each set (INPUT)
C     NSLPTL -- total number of slip systems in all the sets (INPUT)
C     NSET   -- number of sets of slip systems (INPUT)
C
C     H      -- current value of self- and latent-hardening moduli
C               (OUTPUT)
C               H(i,i) -- self-hardening modulus of the ith slip system
C                         (no sum over i)
C               H(i,j) -- latent-hardening molulus of the ith slip
C                         system due to a slip in the jth slip system
C                         (i not equal j)
C
C     PROP   -- material constants characterizing the self- and latent-
C               hardening law (INPUT)
C
C               For the HYPER SECANT hardening law
C               PROP(1,i) -- initial hardening modulus H0 in the ith
C                            set of slip systems
C               PROP(2,i) -- saturation stress TAUs in the ith set of 
C                            slip systems
C               PROP(3,i) -- initial critical resolved shear stress
C                            TAU0 in the ith set of slip systems
C               PROP(9,i) -- ratio of latent to self-hardening Q in the
C                            ith set of slip systems
C               PROP(10,i)-- ratio of latent-hardening from other sets
C                            of slip systems to self-hardening in the
C                            ith set of slip systems Q1
C
C               For Bassani's hardening law
C               PROP(1,i) -- initial hardening modulus H0 in the ith
C                            set of slip systems
C               PROP(2,i) -- stage I stress TAUI in the ith set of 
C                            slip systems (or the breakthrough stress
C                            where large plastic flow initiates)
C               PROP(3,i) -- initial critical resolved shear stress
C                            TAU0 in the ith set of slip systems
C               PROP(4,i) -- hardening modulus during easy glide Hs in
C                            the ith set of slip systems
C               PROP(5,i) -- amount of slip Gamma0 after which a given
C                            interaction between slip systems in the
C                            ith set reaches peak strength
C               PROP(6,i) -- amount of slip Gamma0 after which a given
C                            interaction between slip systems in the
C                            ith set and jth set (i not equal j)
C                            reaches peak strength
C               PROP(7,i) -- representing the magnitude of the strength
C                            of interaction in the ith set of slip
C                            system
C               PROP(8,i) -- representing the magnitude of the strength
C                            of interaction between the ith set and jth
C                            set of system
C               PROP(9,i) -- ratio of latent to self-hardening Q in the
C                            ith set of slip systems
C               PROP(10,i)-- ratio of latent-hardening from other sets
C                            of slip systems to self-hardening in the
C                            ith set of slip systems Q1
C
C     ND     -- leading dimension of arrays defined in subroutine UMAT
C               (INPUT)

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
CC      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL HSELF, HLATNT
CFIXA
      DIMENSION GAMMA(NSLPTL), TAUSLP(NSLPTL), GMSLTL(NSLPTL),
     2          GSLIP(NSLPTL), NSLIP(NSET), PROP(16,NSET),
     3          H(ND,NSLPTL)
CFIXB

      CHECK=0.
      DO I=1,NSET
         DO J=4,8
            CHECK=CHECK+ABS(PROP(J,I))
         END DO
      END DO

C-----  CHECK=0   --  HYPER SECANT hardening law
C       otherwise --  Bassani's hardening law

      ISELF=0
      DO I=1,NSET
         ISET=I
         DO J=1,NSLIP(I)
            ISELF=ISELF+1

            DO LATENT=1,NSLPTL
               IF (LATENT.EQ.ISELF) THEN
CFIXA
                  H(LATENT,ISELF)=HSELF(GAMMA,GMSLTL,GAMTOL,NSLPTL,
     2                                  NSET,NSLIP,PROP(1,I),CHECK,
     3                                  ISELF,ISET)
CFIXB
               ELSE
CFIXA
                  H(LATENT,ISELF)=HLATNT(GAMMA,GMSLTL,GAMTOL,NSLPTL,
     2                                   NSET,NSLIP,PROP(1,I),CHECK,
     3                                   ISELF,ISET,LATENT)
CFIXB

               END IF
            END DO

         END DO
      END DO

      RETURN
      END

C-----------------------------------

C-----  Use single precision on cray
CFIXA
           FUNCTION HSELF(GAMMA,GMSLTL,GAMTOL,NSLPTL,NSET,
     2                           NSLIP,PROP,CHECK,ISELF,ISET)
CFIXB

C-----     User-supplied self-hardening function in a slip system

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
CFIXA
           DIMENSION GAMMA(NSLPTL), NSLIP(NSET), PROP(16),
     2               GMSLTL(NSLPTL)
CFIXB

           IF (CHECK.EQ.0.) THEN

C-----  HYPER SECANT hardening law by Asaro, Pierce et al
              TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              HSELF=PROP(1)*TERM2**2

           ELSE

C-----  Bassani's hardening law
CFIXA
              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)

              ID=0
              G=1.
              DO I=1,NSET
                 IF (I.EQ.ISET) THEN
                    GAMMA0=PROP(5)
                    FAB=PROP(7)
                 ELSE
                    GAMMA0=PROP(6)
                    FAB=PROP(8)
                 END IF

                 DO J=1,NSLIP(I)
                    ID=ID+1
                    IF (ID.NE.ISELF) THEN
CFIXA
               G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
            END IF

                 END DO
              END DO

              HSELF=F*G

           END IF

           RETURN
           END

C-----------------------------------

C-----  Use single precision on cray
CFIXA
           FUNCTION HLATNT(GAMMA,GMSLTL,GAMTOL,NSLPTL,NSET,
     2                            NSLIP,PROP,CHECK,ISELF,ISET,LATENT)
CFIXB

C-----     User-supplied latent-hardening function

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CC           IMPLICIT REAL*8 (A-H,O-Z)
CFIXA
           DIMENSION GAMMA(NSLPTL), NSLIP(NSET), PROP(16),
     2               GMSLTL(NSLPTL)
CFIXB

           ILOWER=0
           IUPPER=NSLIP(1)
           IF (ISET.GT.1) THEN
              DO K=2,ISET
                 ILOWER=ILOWER+NSLIP(K-1)
                 IUPPER=IUPPER+NSLIP(K)
              END DO
           END IF

           IF (LATENT.GT.ILOWER.AND.LATENT.LE.IUPPER) THEN
              Q=PROP(9)
           ELSE
              Q=PROP(10)
           END IF

           IF (CHECK.EQ.0.) THEN

C-----  HYPER SECANT hardening law by Asaro, Pierce et al
              TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              HLATNT=PROP(1)*TERM2**2*Q

           ELSE

C-----  Bassani's hardening law
CFIXA
              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)

              ID=0
              G=1.
              DO I=1,NSET
                 IF (I.EQ.ISET) THEN
                    GAMMA0=PROP(5)
                    FAB=PROP(7)
                 ELSE
                    GAMMA0=PROP(6)
                    FAB=PROP(8)
                 END IF

                 DO J=1,NSLIP(I)
                    ID=ID+1
                    IF (ID.NE.ISELF) THEN
CFIXA
               G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
            END IF

                 END DO
              END DO

              HLATNT=F*G*Q

           END IF

           RETURN
           END

C----------------------------------------------------------------------

CFIXA
      SUBROUTINE ITERATION (GAMMA, TAUSLP, GSLIP, GMSLTL, GAMTOL,
     2                      NSLPTL, NSET, NSLIP, ND, PROP, DGAMOD,
     3                      DHDGDG)
CFIXB

C-----  This subroutine generates arrays for the Newton-Rhapson
C     iteration method.

C-----  Users who want to use their own self- and latent-hardening law
C     may change the function subprograms DHSELF (self hardening) and
C     DHLATN (latent hardening).  The parameters characterizing these
C     hardening laws are passed into DHSELF and DHLATN through array
C     PROP.

C-----  Function subprograms:
C
C       DHSELF -- User-supplied function of the derivative of self-
C                 hardening moduli
C
C       DHLATN -- User-supplied function of the derivative of latent-
C                 hardening moduli

C-----  Variables:
C
C     GAMMA  -- shear strain in all slip systems at the start of time
C               step  (INPUT)
C     TAUSLP -- resolved shear stress in all slip systems (INPUT)
C     GSLIP  -- current strength (INPUT)
CFIX  GMSLTL -- total cumulative shear strains on each individual slip system
CFIX            (INPUT)
C     GAMTOL -- total cumulative shear strains over all slip systems
C               (INPUT)
C     NSLPTL -- total number of slip systems in all the sets (INPUT)
C     NSET   -- number of sets of slip systems (INPUT)
C     NSLIP  -- number of slip systems in each set (INPUT)
C     ND     -- leading dimension of arrays defined in subroutine UMAT
C               (INPUT)
C
C     PROP   -- material constants characterizing the self- and latent-
C               hardening law (INPUT)
C
C               For the HYPER SECANT hardening law
C               PROP(1,i) -- initial hardening modulus H0 in the ith
C                            set of slip systems
C               PROP(2,i) -- saturation stress TAUs in the ith set of 
C                            slip systems
C               PROP(3,i) -- initial critical resolved shear stress
C                            TAU0 in the ith set of slip systems
C               PROP(9,i) -- ratio of latent to self-hardening Q in the
C                            ith set of slip systems
C               PROP(10,i)-- ratio of latent-hardening from other sets
C                            of slip systems to self-hardening in the
C                            ith set of slip systems Q1
C
C               For Bassani's hardening law
C               PROP(1,i) -- initial hardening modulus H0 in the ith
C                            set of slip systems
C               PROP(2,i) -- stage I stress TAUI in the ith set of 
C                            slip systems (or the breakthrough stress
C                            where large plastic flow initiates)
C               PROP(3,i) -- initial critical resolved shear stress
C                            TAU0 in the ith set of slip systems
C               PROP(4,i) -- hardening modulus during easy glide Hs in
C                            the ith set of slip systems
C               PROP(5,i) -- amount of slip Gamma0 after which a given
C                            interaction between slip systems in the
C                            ith set reaches peak strength
C               PROP(6,i) -- amount of slip Gamma0 after which a given
C                            interaction between slip systems in the
C                            ith set and jth set (i not equal j)
C                            reaches peak strength
C               PROP(7,i) -- representing the magnitude of the strength
C                            of interaction in the ith set of slip
C                            system
C               PROP(8,i) -- representing the magnitude of the strength
C                            of interaction between the ith set and jth
C                            set of system
C               PROP(9,i) -- ratio of latent to self-hardening Q in the
C                            ith set of slip systems
C               PROP(10,i)-- ratio of latent-hardening from other sets
C                            of slip systems to self-hardening in the
C                            ith set of slip systems Q1
C
C-----  Arrays for iteration:
C
C       DGAMOD (INPUT)
C
C       DHDGDG (OUTPUT)
C

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
CC      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL DHSELF, DHLATN
CFIXA
      DIMENSION GAMMA(NSLPTL), TAUSLP(NSLPTL), GMSLTL(NSLPTL),
     2          GSLIP(NSLPTL), NSLIP(NSET), PROP(16,NSET),
     3          DGAMOD(NSLPTL), DHDGDG(ND,NSLPTL)
CFIXB

      CHECK=0.
      DO I=1,NSET
         DO J=4,8
            CHECK=CHECK+ABS(PROP(J,I))
         END DO
      END DO

C-----  CHECK=0   --  HYPER SECANT hardening law
C       otherwise --  Bassani's hardening law

      ISELF=0
      DO I=1,NSET
         ISET=I
         DO J=1,NSLIP(I)
            ISELF=ISELF+1

            DO KDERIV=1,NSLPTL
               DHDGDG(ISELF,KDERIV)=0.

               DO LATENT=1,NSLPTL
                  IF (LATENT.EQ.ISELF) THEN
CFIXA
                     DHDG=DHSELF(GAMMA,GMSLTL,GAMTOL,NSLPTL,NSET,
     2                           NSLIP,PROP(1,I),CHECK,ISELF,ISET,
     3                           KDERIV)
CFIXB
                  ELSE
CFIXA
                     DHDG=DHLATN(GAMMA,GMSLTL,GAMTOL,NSLPTL,NSET,
     2                           NSLIP,PROP(1,I),CHECK,ISELF,ISET,
     3                           LATENT,KDERIV)
CFIXB
                  END IF

                  DHDGDG(ISELF,KDERIV)=DHDGDG(ISELF,KDERIV)+
     2                                 DHDG*ABS(DGAMOD(LATENT))
               END DO

            END DO
         END DO
      END DO

      RETURN
      END

C-----------------------------------

C-----  Use single precision on cray
CFIXA
           FUNCTION DHSELF(GAMMA,GMSLTL,GAMTOL,NSLPTL,NSET,
     2                            NSLIP,PROP,CHECK,ISELF,ISET,
     3                            KDERIV)
CFIXB

C-----  User-supplied function of the derivative of self-hardening
C     moduli

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CFIXA
           DIMENSION GAMMA(NSLPTL), GMSLTL(NSLPTL),
     2               NSLIP(NSET), PROP(16)
CFIXB

           IF (CHECK.EQ.0.) THEN

C-----  HYPER SECANT hardening law by Asaro, Pierce et al
              TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              TERM3=PROP(1)/(PROP(2)-PROP(3))*SIGN(1.0,GAMMA(KDERIV))
              DHSELF=-2.*PROP(1)*TERM2**2*TANH(TERM1)*TERM3

           ELSE

C-----  Bassani's hardening law
CFIXA
              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              TERM3=(PROP(1)-PROP(4))/(PROP(2)-PROP(3))

              IF (KDERIV.EQ.ISELF) THEN
                 F=-2.*(PROP(1)-PROP(4))*TERM2**2*TANH(TERM1)*TERM3
                 ID=0
                 G=1.
                 DO I=1,NSET
                    IF (I.EQ.ISET) THEN
                       GAMMA0=PROP(5)
                       FAB=PROP(7)
                    ELSE
                       GAMMA0=PROP(6)
                       FAB=PROP(8)
                    END IF

                    DO J=1,NSLIP(I)
                       ID=ID+1
CFIXA
                       IF (ID.NE.ISELF) G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
                    END DO
                 END DO

              ELSE
                 F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)
                 ILOWER=0
                 IUPPER=NSLIP(1)
                 IF (ISET.GT.1) THEN
                    DO K=2,ISET
                       ILOWER=ILOWER+NSLIP(K-1)
                       IUPPER=IUPPER+NSLIP(K)
                    END DO
                 END IF

                 IF (KDERIV.GT.ILOWER.AND.KDERIV.LE.IUPPER) THEN
                    GAMMA0=PROP(5)
                    FAB=PROP(7)
                 ELSE
                    GAMMA0=PROP(6)
                    FAB=PROP(8)
                 END IF

CFIXA
                 TERM4=GMSLTL(KDERIV)/GAMMA0
CFIXB
                 TERM5=2.*EXP(-TERM4)/(1.+EXP(-2.*TERM4))
                 G=FAB/GAMMA0*TERM5**2

              END IF

              DHSELF=F*G

           END IF

           RETURN
           END

C-----------------------------------

C-----  Use single precision on cray
CFIXA
           FUNCTION DHLATN(GAMMA,GMSLTL,GAMTOL,NSLPTL,NSET,
     2                            NSLIP,PROP,CHECK,ISELF,ISET,LATENT,
     3                            KDERIV)
CFIXB

C-----  User-supplied function of the derivative of latent-hardening
C     moduli

C-----  Use single precision on cray
C
         INCLUDE 'VABA_PARAM.INC'
CFIXA
           DIMENSION GAMMA(NSLPTL), GMSLTL(NSLPTL), NSLIP(NSET),
     2               PROP(16)
CFIXB

           ILOWER=0
           IUPPER=NSLIP(1)
           IF (ISET.GT.1) THEN
              DO K=2,ISET
                 ILOWER=ILOWER+NSLIP(K-1)
                 IUPPER=IUPPER+NSLIP(K)
              END DO
           END IF

           IF (LATENT.GT.ILOWER.AND.LATENT.LE.IUPPER) THEN
              Q=PROP(9)
           ELSE
              Q=PROP(10)
           END IF

           IF (CHECK.EQ.0.) THEN

C-----  HYPER SECANT hardening law by Asaro, Pierce et al
              TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              TERM3=PROP(1)/(PROP(2)-PROP(3))*SIGN(1.0,GAMMA(KDERIV))
              DHLATN=-2.*PROP(1)*TERM2**2*TANH(TERM1)*TERM3*Q

           ELSE

C-----  Bassani's hardening law
CFIXA
              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
              TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
              TERM3=(PROP(1)-PROP(4))/(PROP(2)-PROP(3))

              IF (KDERIV.EQ.ISELF) THEN
                 F=-2.*(PROP(1)-PROP(4))*TERM2**2*TANH(TERM1)*TERM3
                 ID=0
                 G=1.
                 DO I=1,NSET
                    IF (I.EQ.ISET) THEN
                       GAMMA0=PROP(5)
                       FAB=PROP(7)
                    ELSE
                       GAMMA0=PROP(6)
                       FAB=PROP(8)
                    END IF

                    DO J=1,NSLIP(I)
                       ID=ID+1
CFIXA
                       IF (ID.NE.ISELF) G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
                    END DO
                 END DO

              ELSE
                 F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)
                 ILOWER=0
                 IUPPER=NSLIP(1)
                 IF (ISET.GT.1) THEN
                    DO K=2,ISET
                       ILOWER=ILOWER+NSLIP(K-1)
                       IUPPER=IUPPER+NSLIP(K)
                    END DO
                 END IF

                 IF (KDERIV.GT.ILOWER.AND.KDERIV.LE.IUPPER) THEN
                    GAMMA0=PROP(5)
                    FAB=PROP(7)
                 ELSE
                    GAMMA0=PROP(6)
                    FAB=PROP(8)
                 END IF
CFIXA
                 TERM4=GMSLTL(KDERIV)/GAMMA0
CFIXB
                 TERM5=2.*EXP(-TERM4)/(1.+EXP(-2.*TERM4))
                 G=FAB/GAMMA0*TERM5**2

              END IF

              DHLATN=F*G*Q

           END IF

           RETURN
           END

C----------------------------------------------------------------------

      SUBROUTINE LUDCMP (A, N, NP, INDX, D)

C-----  LU decomposition

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
      PARAMETER (NMAX=200, TINY=1.0E-20)
      DIMENSION A(NP,NP), INDX(N), VV(NMAX)

      D=1.
      DO I=1,N
         AAMAX=0.

         DO J=1,N
            IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
         END DO

         IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
         VV(I)=1./AAMAX
      END DO

      DO J=1,N
         DO I=1,J-1
            SUM=A(I,J)

            DO K=1,I-1
               SUM=SUM-A(I,K)*A(K,J)
            END DO

            A(I,J)=SUM
         END DO
         AAMAX=0.

         DO I=J,N
            SUM=A(I,J)

            DO K=1,J-1
               SUM=SUM-A(I,K)*A(K,J)
            END DO

            A(I,J)=SUM
            DUM=VV(I)*ABS(SUM)
            IF (DUM.GE.AAMAX) THEN
               IMAX=I
               AAMAX=DUM
            END IF
         END DO

         IF (J.NE.IMAX) THEN
            DO K=1,N
               DUM=A(IMAX,K)
               A(IMAX,K)=A(J,K)
               A(J,K)=DUM
            END DO

            D=-D
            VV(IMAX)=VV(J)
         END IF

         INDX(J)=IMAX
         IF (A(J,J).EQ.0.) A(J,J)=TINY
         IF (J.NE.N) THEN
            DUM=1./A(J,J)
            DO I=J+1,N
               A(I,J)=A(I,J)*DUM
            END DO
         END IF

      END DO

      RETURN
      END

C----------------------------------------------------------------------

      SUBROUTINE LUBKSB (A, N, NP, INDX, B)

C-----  Linear equation solver based on LU decomposition

C-----  Use single precision on cray
C
    INCLUDE 'VABA_PARAM.INC'
      DIMENSION A(NP,NP), INDX(N), B(N)

      II=0
      DO I=1,N
         LL=INDX(I)
         SUM=B(LL)
         B(LL)=B(I)

         IF (II.NE.0) THEN
            DO J=II,I-1
               SUM=SUM-A(I,J)*B(J)
            END DO
         ELSE IF (SUM.NE.0.) THEN
            II=I
         END IF

         B(I)=SUM
      END DO

      DO I=N,1,-1
         SUM=B(I)

         IF (I.LT.N) THEN
            DO J=I+1,N
               SUM=SUM-A(I,J)*B(J)
            END DO
         END IF

         B(I)=SUM/A(I,I)
      END DO

      RETURN
      END

Pages

Recent comments

More comments

Syndicate

Subscribe to Syndicate