FEAP User Forum
FEAP => Programming => Topic started by: Jfrisch on May 29, 2017, 04:07:03 AM
-
Hello,
I've used the FEAPpv 3.1 a while and I am trying to migrate to the newly purchased FEAP 8.4. After solving some troubles I've managed to compile FEAP 8.4. Now I'm trying to continue debugging my user element (somewehere in my code is an error that cause a divergence of my calculations) but unfortunately I got stuck in an infinite loop of calling the ISW-task 26. While using FEAPpv this ISW-task was never called, so I do not have any experiences in using it.
I've checked the way the built-in solid2d- and solid3d-element handles this task and copied the way (calling the pcorner-subroutine). I've also checked that FEAP 8.4 is executing the code-snippet, where my edits take part by doing an output to stdout (both, the code-snippet in my element and the external subroutine pcorner are called regularly). But, as I mentioned before, unfortunately FEAP never goes past this ISW-task 26 (started a calculation earlier this day and it ran the ISW-task 26 for 30 minutes till I manually aborted it).
Anyone had some analogous problems and can tell me what I need to check or how I probably am able to solve this issue?
Greetings,
Jan Frischmann
-
Hard to say why the infinite loop, but if the code runs using the built-in elements it must relate to your user element.
Can you describe your operating system, compiler used and version, and any other information how you build the program that may give us a clue.
Feap 8.4 does require your user element to have an isw=26 part. It can be dummy that is,
if(isw.eq.26) return
If your element uses the structure
go to(1,2,3, .....) isw
then add a return after this for any isw larger than the list
go to (1,2,3, ...,8) isw
return
Somehow, your element must create a recursion that causes the infinite loop.
The call for isw=26 comes from the routine "pextndc.f" in the ./program directory. You can see if there is something happening after the return (or comment the call toe elmlib to see what happens).
Please let us know what happens if you find the difficulty. Otherwise
if you can post your element we can try to see what is happening.
-
Dear Prof. Taylor,
after declaring all the variables and arrays I do need during my element the following code-snipet is coded by me:
! Store all relevant values for later calls
save
! Set all element-temperatures to zero with double-precision
tl = 0.0d0
! ===========================================================================================
! Set flags for analytical solution and debugging mode
! ===========================================================================================
! set debugging-mode for log-file
debug = .true.
! set output of ISW-task on standard-out
ISWdebug = .true.
if( ISWdebug ) then
write(*,*) ''
write(*,*) 'Element: ',n,' isw: ', isw
endif
! list which jump-mark corresponds to which isw-task (first entry corresponds to isw=1, etc.)
go to ( 1, 2, 3, 4, 5, 3 ), isw
5 return
As far as I know the save-statement should only store all the values assigned to the previously declared variables, after that I am initializing vanishing nodal temperatures (I want to code a purely mechanical element), after that I am initializing to boolean flags and I have coded an output on stdout which tells me which element-number and which ISW-task currently is processed. The main part concerning my problem are the last three lines, in which I am defining, that FEAP should jump to the jumpmarks assigned to the ISW-tasks (if ISW equals 1 it should jump to the first jumpmark '1', if ISW equals 2 FEAP should jump to the jumpmark '2', and so on). If I understand your message correctly FEAP should jump to the jumpmark 3 if ISW equals '6' (the last defined jumpmark in my list) and for every value greather than 6 FEAP should proceed with the code and reach the return-statement after jumpmark 5. So at ISW 26 it should reach this return-statement and exit the code. But obviously FEAP doesn't reassign the ISW-flag/-value to an other value than 26 and so the next time my code is called the previously described behaviour occurs another time and this without ever modifying 'ISW'. So it seems to me, that FEAP needs something to be done at ISW 26, but defining a jumpmark for this ISW (defining jumpmark 5 for every ISW in between 7 and 25) with a code at jumpmark 26 nalogous to the code-snippet of solid2d.f (calling 'pcorner2d()' ) doesn't change anything the same infinite loop occurs.
I have compiled the FEAPpv 3.1 with exact the same code (opened my code for FEAP 8.4 and saved it in the FEAPpv 3.1-user-subfolder) and the code works fine never trying to call ISW 26 (Starting with ISW 1-> 14 -> 12 -> 3). Thats why I am confused of the behaviour of FEAP 8.4, because I have coded the element with the user and programmer manual of FEAP 8.4 so if it works for FEAPpv 3.1 it should work also for FEAP 8.4 (but not necessarily vice versa).
-
Why don't you start the code in gdb and then figure out where it is in the stack when you are stuck in the 'infinite loop'?
Also the behavior of a blanket save statement is to preserve the values of ALL local variables upon exit from a subroutine.
-
Remember that feappv3.1 never has a call to the element with isw = 26, only feap.
Do you use "implicit none" in your element to ensure type casting is correct?
Does feap run correctly with the built in elements?
-
"implicite none" is actually the first line of code of my element and I've tried to run feap with my input file calling for the built-in element solid2d and it exits normally (throwing an IEEE_DIVIDE_BY_ZERO-Exception that isn't thrown by feappv3.1).
-
The IEEE warning can be ignored (it is related to how version 8.4 does NaN and Inf checks).
The main issue is to figure out where you element is hanging. Perhaps some print statements at the top of
the element could help figure out what is going on. Is your element being called in an infinite loop?
-
I have inserted another output to stdout and edited my isw-jump to:
write (*,*) 'Checkpoint 1'
! list which jump-mark corresponds to which isw-task (first entry corresponds to isw=1, etc.)
go to ( 1, 2, 3, 4, 5, 3 ), isw
write (*,*) 'Checkpoint 2'
5 return
write (*,*) 'Checkpoint 3' The first and second Checkpoint are being reached infinitely often but Checkpoint 3 is never reached. So Calling ISW 26 results in executing the return-statement after the go-to-statement but as I said before ISW isn't modified and so my element is being called infinitely often.
-
I've tried the same input file with another user element of mine and the same problem occurs. I've attached my input file and the standard elastic user element with which the loop occurs. Still can't see any problems in my code yet.
-
I do not know what is in your "use math" files, but could they conflict with feap?
I have not compiled because of "use" but make some comments in the file attached (name change to avoid conflict with your file).
Also I would remove all the file I/O to unit 78 until you resolve what is causing the problem.
IF you e-mail the module for "math" we can check further.
-
Just look for RLT in the file to see my comments
-
I compiled your routine with ver84, and it seems to run without infinite looping. The call sequence printed to the screen is isw = 1, 16 x 26, 14, 6, 12, 3, 4. Followed by normal termination. Certainly no infinite loops on isw.eq.26. Note I commented out all the 'use math' lines since you did not provide the module along with all the lines with calls to det( ), inverse( ), and trace( ) which appear to be in that module.
Perhaps you do not have a clean copy of version 8.4?
Note, I use the Intel Fortran Compiler.
-
Oh sorry, was some month ago I used the attached element for the last time, so I forgot about the module I don't use anymore in the new element. The math-module is attached to this post and is only used for distributing some mathematical functions.
-
I compiled your routine with ver84, and it seems to run without infinite looping. The call sequence printed to the screen is isw = 1, 16 x 26, 14, 6, 12, 3, 4. Followed by normal termination. Certainly no infinite loops on isw.eq.26. Note I commented out all the 'use math' lines since you did not provide the module along with all the lines with calls to det( ), inverse( ), and trace( ) which appear to be in that module.
Perhaps you do not have a clean copy of version 8.4?
Note, I use the Intel Fortran Compiler.
I have downloaded the source some days ago because we purchased feap only two weeks ago. I am using the GFortran-Compiler and will try the new release FEAP8.5 with a new and clean installation, checking wether this works.
-
Still resulting in an infinite-loop of calling ISW 26 when compiling the newly released FEAP 8.5. To use my code, which is written in Fortran90, I have adjusted the makefile.in and the makefile in the user-subdirectory (see attachement). No errors printed on stdout and feap tells me that my element (elmt03.f90) is compiled. The benchmark inputfile I am using is also attached to this post. Are there any adjustments I forget about or any issues probably caused by my adjustments?
@Prof. Taylor: You commented in my example-routine, that I do not need to use the ud-array. How can I use the d-array with my user-specific material parameters (how to adjust input-file and in which entries of the d-array are they stored)?
-
Let us avoid calling with isw=26 to see if the error occurs somewhere else too.
In the ./program directory edit pextnd.f to return without doing anything so before the "if(tspine) return" just put "return"
The error does not occur using your element on our systems here. We commented the "use math" and functions so are not doing any real computation.
As for getting data in "ud" it is best to call "modlfd" and "inmate" to use the user material models. Otherwise you will have to follow through all the calls to see how they are obtained from the "d" array that is passed to your element. Start by looking in "modlfd" and find where the user models are done. You should find what the offset is on the "d" to use what we call "ud".
-
I have compiled your routines (with your math.f90) using gfortran in version 8.4. The program runs fine. The output file ends with
*End of Solution Execution* t= 0.00 0.01Here is the isw printed output.
isw: 1
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 26
isw: 14
isw: 6
isw: 12
isw: 3
isw: 4
-
Let us avoid calling with isw=26 to see if the error occurs somewhere else too.
In the ./program directory edit pextnd.f to return without doing anything so before the "if(tspine) return" just put "return"
The error does not occur using your element on our systems here. We commented the "use math" and functions so are not doing any real computation.
I have done that but unluckily there is another strange problem after editing this. I copied the code of my element from feappv3.1-user-subdirectory to the user-subdirectory of FEAP8.5 and compiled FEAP without getting any error-messages. Although code is running for feappv it isn't for FEAP and at the start of the calculation during ISW 6 this error occurs:
Time: 0.0000000000000000 Element: 1 isw: 6
*** Error in `main/feap': free(): invalid pointer: 0x000000000149bac0 ***
======= Backtrace: =========
/lib/x86_64-linux-gnu/libc.so.6(+0x777e5)[0x7f8bd3e687e5]
/lib/x86_64-linux-gnu/libc.so.6(+0x7fe0a)[0x7f8bd3e70e0a]
/lib/x86_64-linux-gnu/libc.so.6(cfree+0x4c)[0x7f8bd3e7498c]
main/feap[0x4bc80d]
main/feap[0x41c61e]
main/feap[0x5462cf]
main/feap[0x40b3a0]
main/feap[0x51b31b]
main/feap[0x40f567]
main/feap[0x4074da]
main/feap[0x4030a8]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf0)[0x7f8bd3e11830]
main/feap[0x4030d9]
======= Memory map: ========
00400000-007ae000 r-xp 00000000 08:03 16781580 /home/jf/FEM/feap85/main/feap
009ad000-009ae000 r--p 003ad000 08:03 16781580 /home/jf/FEM/feap85/main/feap
009ae000-009b2000 rw-p 003ae000 08:03 16781580 /home/jf/FEM/feap85/main/feap
009b2000-00e83000 rw-p 00000000 00:00 0
01465000-014c7000 rw-p 00000000 00:00 0 [heap]
7f8bcc000000-7f8bcc021000 rw-p 00000000 00:00 0
7f8bcc021000-7f8bd0000000 ---p 00000000 00:00 0
7f8bd3113000-7f8bd3118000 r-xp 00000000 08:01 7743199 /usr/lib/x86_64-linux-gnu/libXdmcp.so.6.0.0
7f8bd3118000-7f8bd3317000 ---p 00005000 08:01 7743199 /usr/lib/x86_64-linux-gnu/libXdmcp.so.6.0.0
7f8bd3317000-7f8bd3318000 r--p 00004000 08:01 7743199 /usr/lib/x86_64-linux-gnu/libXdmcp.so.6.0.0
7f8bd3318000-7f8bd3319000 rw-p 00005000 08:01 7743199 /usr/lib/x86_64-linux-gnu/libXdmcp.so.6.0.0
7f8bd3319000-7f8bd331b000 r-xp 00000000 08:01 7743188 /usr/lib/x86_64-linux-gnu/libXau.so.6.0.0
7f8bd331b000-7f8bd351b000 ---p 00002000 08:01 7743188 /usr/lib/x86_64-linux-gnu/libXau.so.6.0.0
7f8bd351b000-7f8bd351c000 r--p 00002000 08:01 7743188 /usr/lib/x86_64-linux-gnu/libXau.so.6.0.0
7f8bd351c000-7f8bd351d000 rw-p 00003000 08:01 7743188 /usr/lib/x86_64-linux-gnu/libXau.so.6.0.0
7f8bd351d000-7f8bd355b000 r-xp 00000000 08:01 7733423 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f8bd355b000-7f8bd375a000 ---p 0003e000 08:01 7733423 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f8bd375a000-7f8bd375b000 r--p 0003d000 08:01 7733423 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f8bd375b000-7f8bd375c000 rw-p 0003e000 08:01 7733423 /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f8bd375c000-7f8bd37ca000 r-xp 00000000 08:01 7753252 /usr/lib/libblas/libblas.so.3.6.0
7f8bd37ca000-7f8bd39c9000 ---p 0006e000 08:01 7753252 /usr/lib/libblas/libblas.so.3.6.0
7f8bd39c9000-7f8bd39ca000 r--p 0006d000 08:01 7753252 /usr/lib/libblas/libblas.so.3.6.0
7f8bd39ca000-7f8bd39cb000 rw-p 0006e000 08:01 7753252 /usr/lib/libblas/libblas.so.3.6.0
7f8bd39cb000-7f8bd39ce000 r-xp 00000000 08:01 5243030 /lib/x86_64-linux-gnu/libdl-2.23.so
7f8bd39ce000-7f8bd3bcd000 ---p 00003000 08:01 5243030 /lib/x86_64-linux-gnu/libdl-2.23.so
7f8bd3bcd000-7f8bd3bce000 r--p 00002000 08:01 5243030 /lib/x86_64-linux-gnu/libdl-2.23.so
7f8bd3bce000-7f8bd3bcf000 rw-p 00003000 08:01 5243030 /lib/x86_64-linux-gnu/libdl-2.23.so
7f8bd3bcf000-7f8bd3bf0000 r-xp 00000000 08:01 7744609 /usr/lib/x86_64-linux-gnu/libxcb.so.1.1.0
7f8bd3bf0000-7f8bd3def000 ---p 00021000 08:01 7744609 /usr/lib/x86_64-linux-gnu/libxcb.so.1.1.0
7f8bd3def000-7f8bd3df0000 r--p 00020000 08:01 7744609 /usr/lib/x86_64-linux-gnu/libxcb.so.1.1.0
7f8bd3df0000-7f8bd3df1000 rw-p 00021000 08:01 7744609 /usr/lib/x86_64-linux-gnu/libxcb.so.1.1.0
7f8bd3df1000-7f8bd3fb0000 r-xp 00000000 08:01 5243034 /lib/x86_64-linux-gnu/libc-2.23.so
7f8bd3fb0000-7f8bd41b0000 ---p 001bf000 08:01 5243034 /lib/x86_64-linux-gnu/libc-2.23.so
7f8bd41b0000-7f8bd41b4000 r--p 001bf000 08:01 5243034 /lib/x86_64-linux-gnu/libc-2.23.so
7f8bd41b4000-7f8bd41b6000 rw-p 001c3000 08:01 5243034 /lib/x86_64-linux-gnu/libc-2.23.so
7f8bd41b6000-7f8bd41ba000 rw-p 00000000 00:00 0
7f8bd41ba000-7f8bd41d0000 r-xp 00000000 08:01 5247454 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8bd41d0000-7f8bd43cf000 ---p 00016000 08:01 5247454 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8bd43cf000-7f8bd43d0000 rw-p 00015000 08:01 5247454 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8bd43d0000-7f8bd44d8000 r-xp 00000000 08:01 5243032 /lib/x86_64-linux-gnu/libm-2.23.so
7f8bd44d8000-7f8bd46d7000 ---p 00108000 08:01 5243032 /lib/x86_64-linux-gnu/libm-2.23.so
7f8bd46d7000-7f8bd46d8000 r--p 00107000 08:01 5243032 /lib/x86_64-linux-gnu/libm-2.23.so
7f8bd46d8000-7f8bd46d9000 rw-p 00108000 08:01 5243032 /lib/x86_64-linux-gnu/libm-2.23.so
7f8bd46d9000-7f8bd4802000 r-xp 00000000 08:01 7744222 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f8bd4802000-7f8bd4a01000 ---p 00129000 08:01 7744222 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f8bd4a01000-7f8bd4a02000 r--p 00128000 08:01 7744222 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f8bd4a02000-7f8bd4a04000 rw-p 00129000 08:01 7744222 /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f8bd4a04000-7f8bd4a5b000 r-xp 00000000 08:01 7743831 /usr/lib/x86_64-linux-gnu/libjpeg.so.8.0.2
7f8bd4a5b000-7f8bd4c5b000 ---p 00057000 08:01 7743831 /usr/lib/x86_64-linux-gnu/libjpeg.so.8.0.2
7f8bd4c5b000-7f8bd4c5c000 r--p 00057000 08:01 7743831 /usr/lib/x86_64-linux-gnu/libjpeg.so.8.0.2
7f8bd4c5c000-7f8bd4c5d000 rw-p 00058000 08:01 7743831 /usr/lib/x86_64-linux-gnu/libjpeg.so.8.0.2
7f8bd4c5d000-7f8bd5251000 r-xp 00000000 08:01 7753255 /usr/lib/lapack/liblapack.so.3.6.0
7f8bd5251000-7f8bd5450000 ---p 005f4000 08:01 7753255 /usr/lib/lapack/liblapack.so.3.6.0
7f8bd5450000-7f8bd5451000 r--p 005f3000 08:01 7753255 /usr/lib/lapack/liblapack.so.3.6.0
7f8bd5451000-7f8bd5455000 rw-p 005f4000 08:01 7753255 /usr/lib/lapack/liblapack.so.3.6.0
7f8bd5455000-7f8bd558a000 r-xp 00000000 08:01 7743186 /usr/lib/x86_64-linux-gnu/libX11.so.6.3.0
7f8bd558a000-7f8bd578a000 ---p 00135000 08:01 7743186 /usr/lib/x86_64-linux-gnu/libX11.so.6.3.0
7f8bd578a000-7f8bd578b000 r--p 00135000 08:01 7743186 /usr/lib/x86_64-linux-gnu/libX11.so.6.3.0
7f8bd578b000-7f8bd578f000 rw-p 00136000 08:01 7743186 /usr/lib/x86_64-linux-gnu/libX11.so.6.3.0
7f8bd578f000-7f8bd57b5000 r-xp 00000000 08:01 5243035 /lib/x86_64-linux-gnu/ld-2.23.so
7f8bd5984000-7f8bd598c000 rw-p 00000000 00:00 0
7f8bd59b1000-7f8bd59b4000 rw-p 00000000 00:00 0
7f8bd59b4000-7f8bd59b5000 r--p 00025000 08:01 5243035 /lib/x86_64-linux-gnu/ld-2.23.so
7f8bd59b5000-7f8bd59b6000 rw-p 00026000 08:01 5243035 /lib/x86_64-linux-gnu/ld-2.23.so
7f8bd59b6000-7f8bd59b7000 rw-p 00000000 00:00 0
7ffdfcb91000-7ffdfcbb2000 rw-p 00000000 00:00 0 [stack]
7ffdfcbe8000-7ffdfcbea000 r--p 00000000 00:00 0 [vvar]
7ffdfcbea000-7ffdfcbec000 r-xp 00000000 00:00 0 [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0 [vsyscall]
Program received signal SIGABRT: Process abort signal.
Backtrace for this error:
#0 0x7F8BD46F4E08
#1 0x7F8BD46F3F90
#2 0x7F8BD3E264AF
#3 0x7F8BD3E26428
#4 0x7F8BD3E28029
#5 0x7F8BD3E687E9
#6 0x7F8BD3E70E09
#7 0x7F8BD3E7498B
#8 0x4BC80C in elmt03_ at elmt03.f90:49 (discriminator 1)
#9 0x41C61D in elmlib_ at elmlib.f:64
#10 0x5462CE in pform_ at pform.f:338
#11 0x40B39F in formfe_ at formfe.f:92
#12 0x51B31A in pmacr2_ at pmacr2.f:599
#13 0x40F566 in pmacr_ at pmacr.f:591
#14 0x4074D9 in pcontr_ at pcontr.f:1375
#15 0x4030A7 in feap at feap85.f:178
#16 0x7F8BD3E1182F
Aborted (core dumped)
The first line is an output generated because of my code and shows which ISW-task FEAP is currently processing. Following the backtrace leads to the 49th line of my element and as you can see in the following code-snippet, line 49 of my code is the allocation of the shp-array. Is there any obvious fault in my code that I can't see? Changing the flexible dimensions "ndm+1" and "nen" to fixed values results in other errors depending on the value I pass there.
subroutine elmt03(d,ul,xl,ix,tl,stiff,residual,ndf,ndm,nst,isw)
! ============================================================================================
! * * F E A P * * A FiniteElement Analysis Program
! .... Copyright (c) 1984-2013: Robert L. Taylor
! All rights reserved
! ============================================================================================
implicit none
include 'eldata.h' ! n_el - current element
! ma - material set of current element
! mct - print counter
! iel - user element number
! nel - number of nodes of current element
include 'cdata.h' ! numnp - number of mesh nodes
! numel - number of mesh elements
! nummat - number of material sets
! nen - maximum number of nodes per element
! neq - number of equations
! ipr - real variable precision
include 'tdata.h' ! ttim - actual time
! dt - time-increment
include 'hdata.h' ! nhi - history-variables 1,2,3
include 'eltran.h' ! ctan-modifiers
! Inputs
real*8 :: d(*) ! FEAP-internal material parameters
real*8 :: ul(ndf,nen,*) ! nodal displacement vector
! (third index is '1' for time n+1 & '2' for increment)
real*8 :: xl(ndm,nen) ! nodal coordinates of elements
real*8 :: tl(nen) ! nodal temperatures
integer :: ix(nen) ! global node number of local nodes -> ID-Array
integer :: ndf ! number of degrees of freedom per node
integer :: ndm ! number of dimensions of mesh/problem
integer :: nst ! number of stresses = ndf * nen
integer :: isw ! control parameter for element routine
! Outputs
real*8 :: stiff(nst,nst) ! stiffness-matrix (FEAP-internal: s)
real*8 :: residual(nst) ! residual force vector (FEAP-internal: p)
! Internal calculation parameters
integer,parameter :: nh=1 ! number of history variables for each time step
integer :: lint ! number of quadrature points
integer :: tdof ! thermal degrees of freedom
real*8 :: coor(ndm) ! temporary array for coordinates of Gauss-point
real*8 :: shp(ndm+1,nen) ! Shape-Functions and their derivatives
real*8 :: gauss(ndm+1,nen) ! coordinates and weights of Gauss-points
real*8 :: ud(10) ! user-defined material parameters
real*8 :: dV ! volume increment in reference configuration
real*8 :: dV_e ! volume of current element (approximation)
real*8 :: jacobian ! Jacobian-Matrix
! counting indices
integer :: i,j,k,l,m,n,pp ! counting indices
integer :: A,B,C,DD ! counting indices for nodes
integer :: P,Q ! multi indices for microdeformation
integer :: G ! counting index for Gauss-Iteration
integer :: count ! counting number for ISW-3 task
! booleans for flags of calculation
logical :: debug ! flag for debugging
logical :: ISWdebug ! flag for output of ISW
logical :: analytical ! flag for analytical material tangent moduli
logical :: errck, pinput
! scalar valued material parameters
real*8 :: mu_1
real*8 :: lambda
real*8 :: c_1
real*8 :: c_3
! kinematical tensors
real*8 :: H(3,3) ! displacement gradient
real*8 :: E(3,3) ! linearized Green's strain tensor
real*8 :: Psi(3,3) ! microdeformation
real*8 :: dPsi(3,3,3) ! gradient of microdeformation
! stress tensors
real*8 :: Sigma(3,3) ! Cauchy stresses
real*8 :: tau(3,3) !
real*8 :: mu(3,3,3) ! microstresses
real*8 :: Sigma_e(3,3) ! stress tensor of element (approximation)
! elasticity tensor
real*8 :: Ecoeff(3,3,3,3)
! tensorial quantities for FEM-stiffness
real*8 :: Chi_1(3,3,3,3) ! coefficient of K_UU
real*8 :: Chi_2(3,3,3,3) ! coefficient of K_UPsi
real*8 :: Chi_3(3,3,3,3) ! coefficient of K_PsiU
real*8 :: Chi_4(3,3,3,3,3,3) ! second coefficient of K_PsiPsi
! mathematical tensors
real*8 :: kron(3,3) ! Kronecker delta
real*8 :: perm(3,3,3) ! permutation tensor
! multi index for stiffness of microdeformation
integer :: MultiIndex(6,2)
character*100 datfile
! Store all relevant values for later calls
save
Are there any major changes in passing the arguments to the shp2d or int2d-subroutines made? Or are there some major changes how these two FEAP-subroutines work with respect to the FEAPpv-version?
-
Sounds like you have included the incorrect include/integerX files, where X=4 or 8.
-
If you compile maintain/memprobe.c and run it, you will get an output that tells you the details of your computer and which include files to pick as well as how to select the value of ipr in main/feap85.f.
-
**** FEAP Memory Prober: Start ****
Your machine's data-types use the following number of bits
char 8, int 32, long 64, long long 64, float 32, double 64, pointers 64
Select the integer8 include files in the top level makefile.in
Set ipr to 2 in main/feapXX.f
**** FEAP Memory Prober: End ****
That is the output of a.out and as you can see in the attached files I have chosen the integer8-include-directory for a 64bit-PC and ipr=2.
-
This looks fine but the free error that you posted is 99% of the time due to the incorrect selection of pointer sizes. Can you try deleting
all your *.o and *.a files and re-building?
find $FEAPHOME8_5 -name "*.o *.a" -exec rm {} \;
If that fails to fix the problem, please e-mail your files to feap@berkeley.edu .
-
I think I have figured out the latest problem and it is the flexible declaration of two arrays in my code (ndm and nen is set during runtime and after replacing these boundaries with fixed values the pointer exception vanished). But now there is another problem that occured after I had down some minor fixes to my feappv code (and testes it with feappv). When I trap excpetions by adding the "-ffpe-trap=overflow,underflow,zero"-compileroption I get a "erroneous arithmetic operation"-exception because there is a division by zero. I have traced the failure by the backtrace on stdout and the division by zero occurs in the inmate.f-subroutine at line 2335:
dd(3,3) = 1.0d0/e3
I have tried to figure out what causes this problem and I have no clue (even using a debugger proceeding step by step I didn't see when actually the value e3 is set to zero). The cod-snippet that is used of my own element (because the exception is thrown after calling ISW=1) is the following:
subroutine elmt03(d,ul,xl,ix,tl,stiff,residual,ndf,ndm,nst,isw)
! ============================================================================================
! * * F E A P * * A FiniteElement Analysis Program
! .... Copyright (c) 1984-2013: Robert L. Taylor
! All rights reserved
! ============================================================================================
implicit none
include 'eldata.h' ! n_el - current element
! ma - material set of current element
! mct - print counter
! iel - user element number
! nel - number of nodes of current element
include 'cdata.h' ! numnp - number of mesh nodes
! numel - number of mesh elements
! nummat - number of material sets
! nen - maximum number of nodes per element
! neq - number of equations
! ipr - real variable precision
include 'tdata.h' ! ttim - actual time
! dt - time-increment
include 'hdata.h' ! nhi - history-variables 1,2,3
include 'eltran.h' ! ctan-modifiers
! Inputs
real*8 :: d(*) ! FEAP-internal material parameters
real*8 :: ul(ndf,nen,*) ! nodal displacement vector
! (third index is '1' for time n+1 & '2' for increment)
real*8 :: xl(ndm,nen) ! nodal coordinates of elements
real*8 :: tl(nen) ! nodal temperatures
integer :: ix(nen) ! global node number of local nodes -> ID-Array
integer :: ndf ! number of degrees of freedom per node
integer :: ndm ! number of dimensions of mesh/problem
integer :: nst ! number of stresses = ndf * nen
integer :: isw ! control parameter for element routine
! Outputs
real*8 :: stiff(nst,nst) ! stiffness-matrix (FEAP-internal: s)
real*8 :: residual(nst) ! residual force vector (FEAP-internal: p)
! Internal calculation parameters
integer,parameter :: nh=1 ! number of history variables for each time step
integer :: lint ! number of quadrature points
integer :: tdof ! thermal degrees of freedom
real*8 :: coor(ndm) ! temporary array for coordinates of Gauss-point
real*8 :: shp(4,4) ! Shape-Functions and their derivatives
real*8 :: gauss(4,4) ! coordinates and weights of Gauss-points
real*8 :: ud(10) ! user-defined material parameters
real*8 :: dV ! volume increment in reference configuration
real*8 :: jacobian ! Jacobian-Matrix
! counting indices
integer :: i,j,k,l,m,n,pp ! counting indices
integer :: A,B ! counting indices for nodes
integer :: P,Q ! multi indices for microdeformation
integer :: G ! counting index for Gauss-Iteration
integer :: count ! counting number for ISW-3 task
! booleans for flags of calculation
logical :: debug ! flag for debugging
logical :: ISWdebug ! flag for output of ISW
logical :: analytical ! flag for analytical material tangent moduli
logical :: errck, pinput
! scalar valued material parameters
real*8 :: mu_1
real*8 :: lambda
real*8 :: c_1
real*8 :: c_3
! kinematical tensors
real*8 :: H(3,3) ! displacement gradient
real*8 :: E(3,3) ! linearized Green's strain tensor
real*8 :: Psi(3,3) ! microdeformation
real*8 :: dPsi(3,3,3) ! gradient of microdeformation
! stress tensors
real*8 :: Sigma(3,3) ! Cauchy stresses
real*8 :: tau(3,3) !
real*8 :: mu(3,3,3) ! microstresses
real*8 :: Sigma_e(3,3) ! stress tensor of element (approximation)
! elasticity tensor
real*8 :: Ecoeff(3,3,3,3)
! tensorial quantities for FEM-stiffness
real*8 :: Chi_1(3,3,3,3) ! coefficient of K_UU
real*8 :: Chi_2(3,3,3,3) ! coefficient of K_UPsi
real*8 :: Chi_3(3,3,3,3) ! coefficient of K_PsiU
real*8 :: Chi_4(3,3,3,3,3,3) ! second coefficient of K_PsiPsi
! mathematical tensors
real*8 :: kron(3,3) ! Kronecker delta
real*8 :: perm(3,3,3) ! permutation tensor
! multi index for stiffness of microdeformation
integer :: MultiIndex(6,2)
! Store all relevant values for later calls
save
! Set all element-temperatures to zero with double-precision
tl = 0.0d0
! ===========================================================================================
! Set flags for analytical solution and debugging mode
! ===========================================================================================
! set debugging-mode for log-file
debug = .false.
! set output of ISW-task on standard-out
ISWdebug = .true.
! set analtical material tangent moduli
analytical = .true.
if( ISWdebug ) then
write(*,*) ''
write(*,*) 'Time:',ttim,'Element:',n_el,' isw:', isw
endif
! list which jump-mark corresponds to which isw-task (first entry corresponds to isw=1, etc.)
go to ( 1, 2, 3, 4, 5, 3 ), isw
5 return
! ============================================================================================
! ISW = 1: Input d(*) parameters
! ============================================================================================
! write mesh-information to log-file
1 open( 78, file='log.txt', form='formatted' )
write(78,2000)
write(78,2001) numnp, numel, nen, ndm, ndf, nummat
close(78)
! read thermal degrees of freedom
tdof = nint(d(19))
! call routine for material parameters
call inmate ( d, tdof, nh, 1 )
! Output:
! d -> set of material parameters
! tdof -> thermal degrees of freedom
! nh -> number of history variables used (1)
!
! Input:
! 1 -> type of element (here: mechanical solid)
I pass a newly initiated d-array (declared by "d(*)"), the thermal degrees of freedom (0), the amount of history variables used (I pass a 1 here but I don't use any history variable so 0 should be ok too) and I declare that I am using a mechanical solid by the last argument "1". I don't see any faults in my coding there especially because feappv isn't rumoring in any manner. So did you change the way the inmate-subroutine operates and I haven't figured out the change or can you see or guess another problem here?
-
This is telling you that you have not properly set your elastic properties.
Have a look at the line you are using to set your elastic moduli; there is something wrong
with it. Also look in inmate.f at the section being used for your material model.
-
I have used the attached input-file (working without problem with feappv) ... I don't see the problem (I use the ud-Array for my elastic material parameters, haven't figured out yet how to use the d-array for user-defined parameters).
-
Why the blank line after ucon?
You should not have to use the 'n1' and 'e1' variables feap should be able to count how many were in each block.
Do not try to run the problem after the tie. Do a graphical check that you got everything in correctly. Check the figure of the mesh, the boundary conditions, the loads if any, etc. Also check the "O" and "L" files for any warnings or errors. After you have checked that the problem looks o.k. then try one solution to see if there are any solver or other warnings, you may also check plots of contours, etcl -- make sure there are no rigid body modes, or other issues. Once you are confident things look o.k. try to converge the first step -- then do the problem.
-
1.) I added the blank line because my colleagues told me to (I used/looked at their input files in first place as a reference). But deleted it now.
2.) I am using the variables "e1" and "n1" for the next block, at which I need to assign a value to the first number of node and element according to the manual. If this is obsolet I don't understand why these parameters are coded?
3.) Both, the "O"- and the "L"-file, are empty, so there are no warnings or errors. If I turn off the ffpe-trap "zero" I ionly get this error-message:
*D4TRI WARNING* Reduced diagonal is zero in 14742 equations.
I would suggest that this warning is according to the exception thrown if I turn the ffpe-trap on. The boundary-conditions and loads look fine in the "O"-file (exactly what I tried to do with my input-file). If I edit my input-file to the following code:
FEAP * * 4-Element Patch Test
0 0 0 2 5 4
! total amount of nodes (calculated automatically)
! total amount of elements (calculated automatically)
! total amount of materials used (calculated automatically)
! 2 dimensions (mesh)
! 5 degrees of freedom (mesh)
! 4 nodes (per element)
PARAmeters
n1 = 1 ! index of master (first) node of element
e1 = 1 ! index of element
xi = 100 ! number of elements generated along radial axis in inner square
yi = 15 ! number of elements generated along spherical axis
l1 = 0.25 ! radius of hole
l2 = 2.5 ! one half of length of inner square
l3 = l1/sqrt(2) ! coordinates of point half-way of arc
l4 = cosd(22.5)*l1 ! x-coordinate of point at one forth of arclength
l5 = sind(22.5)*l1 ! y-coordinate of point at one forth of arclength
l6 = (l2+2*l1)/3 ! intermdiate point along horizontal/vertical axes of inner square
l7 = (l2+2*l3)/3 ! intermediate point along the diagonal axes of inner square
e = 210000 ! Young's modulus
v = 0.3 ! Poission's ratio
c1 = 1
c3 = 5
u2 = 0.05 ! boundary condition for upper edge
MATErial ! specify material used
USER 3
UCON mat1
e v c1 c3
! blank termination record
BLOCk ! north northeast
CARTesian xi yi n1 e1 1
1 l3 l3
2 l2 l2
3 0 l2
4 0 l1
5 l7 l7
7 0 l6
8 l5 l4
! termination record
PARAmeters ! adjust numbering of nodes and elements
n1 = n1 + (xi+1)*(yi+1)
e1 = e1 + xi*yi
BLOCk ! east northeast
CARTesian xi yi n1 e1 1
1 l1 0
2 l2 0
3 l2 l2
4 l3 l3
5 l6 0
7 l7 l7
8 l4 l5
! termination record
EBOUndary ! set boundary conditions for symmetric mesh
2 l2 0 1 1 1 1 ! u_y and psi vanishing for upper edge
2 0 0 1 1 1 1 ! vanishing u_y and psi for lower edge
1 l2 0 0 1 1 1 ! vanishing psi for right edge
1 0 1 0 1 1 1 ! vanishing u_x and psi for left edge
CBOUndary ! set additional boundary conditions for lower right corner
node 0 l2 1 1 1 1 1
EDISplacement ! set boundary conditions (amount)
2 l2 0 u2 0 0 0
! termination record
END
TIE ! tie all nodes with same coordinates (within tolerance) together
inte
end
! batch
! prop,,3
! dt,,0.1
! tol,,1e-10
! loop,,10
! time
! loop,,5
! UTANg,,1
! next
! disp all
! stre all
! next
! end
! 2,2
! 0.,0. 1.,1.
! STOP
I get exactly the same error-message respectively exception thrown (output of stdout) without starting the FEAP-GUI to check mesh with FEAP:
F I N I T E E L E M E N T A N A L Y S I S P R O G R A M
FEAP (C) Regents of the University of California
All Rights Reserved.
VERSION: Release 8.5.2a
DATE: 26 May 2017
Files are set as: Status Filename
Input (read ) : Exists Iloch
Output (write) : Exists Oloch
Restart (read ) : New Rloch
Restart (write) : New Rloch
Plots (write) : New Ploch
Caution, existing write files will be overwritten.
Are filenames correct?( y or n; r = redefine all, s = stop) : y
R U N N I N G F E A P P R O B L E M N O W
--> Please report errors by e-mail to:
feap@berkeley.edu
isw: -1
Time: 0.0000000000000000 Element: 0 isw: -1
Time: 0.0000000000000000 Element: 0 isw: 1
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 0x7F9092C33E08
#1 0x7F9092C32F90
#2 0x7F90923654AF
#3 0x45FB79 in inmate_ at inmate.f:2335
#4 0x4BC99E in elmt03_ at elmt03.f90:142
#5 0x41C67D in elmlib_ at elmlib.f:64
#6 0x6C9F29 in pmatin_ at pmatin.f:280
#7 0x5BAB8A in pmesh_ at pmesh.f:464
#8 0x4F5B86 in pnewprob_ at pnewprob.f:740
#9 0x40442F in pcontr_ at pcontr.f:1009
#10 0x403101 in feap at feap85.f:178
#11 0x7F909235082F
Floating point exception (core dumped)
I have checked the mesh in feappv before starting calculations with it and the mesh generated looks fine too. There are now degenerated elements visible and as I mentioned earlier the boundary conditions are set correct.
-
1. The way to conceptualize the need for blank lines is that that is how many of FEAP's internal command decide when the input they are supposed to read ends. Thus blank lines are like semi-colons (;) in C programs or </xxxx> tags in markup languages.
2. You do not need n1 e1, FEAP will auto count for you if you leave off the n1 e1 from the CARTesian line; just give CARTesian x1 y1. FEAP will count up the elments and nodes for you and set the material to the default (first) material.
3. The datri warning is telling you that you have singular equations. Someplace your input is not working. You should print out your material properties to make sure they are getting read correctly; check the O file too.
-
1. The way to conceptualize the need for blank lines is that that is how many of FEAP's internal command decide when the input they are supposed to read ends. Thus blank lines are like semi-colons (;) in C programs or </xxxx> tags in markup languages.
2. You do not need n1 e1, FEAP will auto count for you if you leave off the n1 e1 from the CARTesian line; just give CARTesian x1 y1. FEAP will count up the elments and nodes for you and set the material to the default (first) material.
3. The datri warning is telling you that you have singular equations. Someplace your input is not working. You should print out your material properties to make sure they are getting read correctly; check the O file too.
1. Ah tanks, that clearifies the intended behaviour for me.
2. That sounds very nice - I will try and do it that way in the future.
3. As I mentioned above the "O"-file is empty - so there are especially no warning or error messages. Except the datri-warning there are no other warnings or errors occuring in the "O"-file while disabling the ffpe-trap=zero compiler-flag (and BC, loads, elements, material properties look fine). If I enable it the above shown exception is thrown to my stdout and I do not know at which line my input-file is incorrect (because it works fine with feappv and I have build it according to the FEAP 8.4-manual so if there are issues like that I would have guessed them the other way: input file working with feap 8.4 or 8.5 and not working with feappv 3.1 but not vice versa). I have also checked the material properties with feappv and they were fine. The exception of FEAP 8.5 is thrown while these material properties are read from input-file (at least I guess so because the exception is thrown by inmate.f).
-
Why don't you compile with -g and then use a debugger to try and trap the error? or at the minimum step through the code to see where things are going bad?
-
I am using the -g option and already tried to trap the error with a debugger, but the erroneous variable "e3" is never used or set (except the declaration) till the line at which the exception is thrown. I have absolutely no clue what's causing the exception.
-
Maybe a better option will be to simply skip using inmate and do your own material property inputs. When you use your own elements there are no requirements to use inmate.