Output internal /history variables in feap user element routine

[ace21]

Dear feap community

I have a user defined element subroutine, from which I want to print out some internal/history variables

during the computation of tangent sttiffness(isw=3) I compute the updated internal variables and assign them to array hr e.g.

do i=1,5
hr(nh2-1+i)=bet(i)
end do
where bet(i) was computed during tangent stiffnes computation phase

however when I try to see the  output of these values  in the output file, it prints zeros for the corresponding variables in the output file
 e.g.
isw.eq.4
write(iow,*)'sig11',hr(nh2-1+1),'sig22',hr(nh2-1+2),....
....

could you kindly suggest what could possibly be the error

If I am not clear enough I will try to explain in more detail

regards

[Prof. R.L. Taylor]

Did you set nh1 under is isw=1?

[ace21]

yes I did set its value

[FEAP_Admin]

To double check the allocation you should issue 'SHOW,DICTionary' as your very first Macro command.  You can then see if the history array has been properly allocated.

[ace21]

Dear Prof Taylor

I issued the command 'show,dict' after running the simulation I got the following stats for history variables

             30   H            49      2       42           2293349 Program
             31   NH1       50      2       42           2293397 Program
             32   NH2       51      2       42           2293445 Program
I guess the array is properly alocated ( I have in total  21 values that are stored in it)
or kindly correct  me if I have interpreted these values wrongly

[Prof. R.L. Taylor]

How many elements are you using?  It appears only 1 -- is this correct?

If you issue the command SHOW H

you will get all the values that are currently stored in the history array.  When you are running an analysis, the values in H are transferred to the arrays NH1 and NH2 for each individual element and passed to the subroutine using 'comblk.h' and 'hdata.h'.

comblk.h has the pointer to the array HR(*) and hdata.h has the value of nh1 and nh2 that define the offset to where the data starts.  THus, your print should be getting what is available in H.  However, to have feap really save the data, that is transfer back to H what you put into hr(nh2+...) you must be in a part of the program where data is saved.  By default this is only when isw = 3 or 6 or 14.  So if you are computing stresses and storing only in isw=4 the data is not transferred from the hr(nh1..), hr(nh2..) or hr(nh3..).  This is explained in the programmer manual.

The reasons for the restrictions is evident.  If you changed the history values when you were performing an output of variables then when you started the next solution things would be corrupted by what you did during a print only operation.

[ace21]

Thanks for the detailed explanation Prof Taylor

 Yes I am using only one element

[ace21]

From your explaination, the values updated/stored in hr can only be recovered later for plot/print if the update is done in   isw=3,6 or14, using this fact I went through the element routine to see where the storage in hr takes place, here is the essential workflow of the element

if(isw.eq.1) then
nh1 = 21 ( assign number of history terms)
elseif(isw.eq.2)
check mesh
elseif(isw.eq. 3 .or. isw.eq. 4 .or. isw.eq. 6 .or. isw.eq. 8 .or. isw.eq.14) then
assign integration points(lint)
if(isw.eq.14) then
 initialize some of history variables
 do i = 1,lint
 hr(nh1+ni + (i-1)*nhv + 3) = 1.0d0/9.0d0/d(uprm+1)
 hr(nh2+ni + (i-1)*nhv + 3) = hr(nh1+ni + (i-1)*nhv + 3)
 enddo ! i
endif !(isw.eq.14)
get values from hr(nh1 + ..) and hr(nh2+..) for some calculations
e.g stress interpolation parameters
do i = 1,ni
 bet(i)   = hr(nh2-1+i)
 write(*,*)'bet(i)',bet(i)    
 betn(i)  = hr(nh1-1+i)
end do
computing strain,jacobi and other stuff
if(isw.eq.3 .or. isw.eq.6 .or. isw.eq.14) then
if(isw.eq.3 .or. isw.eq.6) then
here the tangent stiffness  and residual is calculated
and history variables at tn+1 are updated (among other variables)e.g. the newly
computed bet(i) (in isw.eq.3 .or. isw.eq.6) is assigned to hr array
 do i = 1,ni
 hr(nh2-1+i) = bet(i)
 end do
endif !( isw = 3 or 6  )
elseif(isw.eq.4 .or. isw.eq. then
get values from hr(nh2+..) for output
and plotting
HERE IS THE PROBLEM: ALL VALUES ARE ZERO in hr
endif !(isw = 4 or
endif !(isw = 3 or 4 or 6 or 8 or 14 or 16)
what I have understood from the code that the storage in hr does take place in (isw = 3 or 6)  and in (isw = 4 or I only retrieve the values in hr(nh2+..) and do not modify them, however they are all zero, even though I had updated them during (isw = 3 or 6) or am I missing some thing in understanding

Thanks for your patience and advice

[Prof. R.L. Taylor]

I assume the bet(i) you output to the screen were non zero?

If you do SHOW H is it all zero?  If not:

Can you post your routine(s) so I can see exactly what is done?

What you describe seem o.k. but somewhere it is being lost before begin put back in the H array storage.

[ace21]

Dear Prof Taylor

I printed out bet(i) it also shows zero values

and SHOW H also shows its empty, I am really confused now whats wrong with the code

Anyway I am attaching the element subroutine, may be you can catch the error

[Prof. R.L. Taylor]

You need to be careful that when you do a computation in the element that all parameters are defined properly.  Do not assume they are still there just because you set something in a previous isw.   For example in isw=14 I see lint but it is still zero probably since you never set it -- so nothing is ever set in the array. 

Try to structure your code better -- indent evenly to see where you are in the logic and computation, etc.

Setting lint should help get things going

[ace21]

Dear Prof Taylor

I printed the value of lint in isw=14  and also 3,6,8 etc, it is set to its proposed value i.e. 4

I agree the code is not well structured, actually it took me a while to understand the work flow since I did not code it myself

Anyways thanks for your valuable advice, I am still new to feap, but this forum is very useful in learning and implementing ur own stuff in feap