Add a New Element in Fortran (FORM_TANG_AND_RISID Block)

[dray]

Hi all,

I follow the instruction on http://opensees.berkeley.edu/wiki/index ... nt_Fortran to build a trussf element. The dll runs successfully. However, I noticed when I added a write command in the block of
ISW = =ISW_ FORM_TANG_AND_RISID as below, it will record the resid and tang twice at each step.

In ISW_FORM_TANG_AND_RESID, all the trials during a non-linear analysis are performed. DO NOT save state variables here.
call c_f_pointer(eleObj%param, theParam, [4]);
call c_f_pointer(eleObj%node, theNodes, [2]);
call c_f_pointer(eleObj%mats, theCMatPtrPtr, [1]);
theCMatPtr = theCMatPtrPtr(1);

A = theParam(1);
L = theParam(2);
cs = theParam(3);
sn = theParam(4);
nd1 = theNodes(1);
nd2 = theNodes(2);

numDOF = 2;
err = OPS_GetNodeDisp(nd1, numDOF, d1);
err = OPS_GetNodeDisp(nd2, numDOF, d2);

tran(1) = -cs;
tran(2) = -sn;
tran(3) = cs;
tran(4) = sn;

dLength = 0.0;
do i = 1,2
dLength = dLength - (d2(i)-d1(i)) * tran(i);
continue

strn(1) = dLength/L;

c i = 0
c i=OPS_InvokeMaterial(eleObj, i, modl, strn, strs, tng, isw)
j=OPS_InvokeMaterialDirectly(theCMatPtr, modl, strn, strs,
+ tng, isw)

force = A*strs(1);
k = A*tng(1)/L;

do i =1,4
resid(i) = force * tran(i);
do j = 1,4
tang(i,j) = k * tran(i)*tran(j);
continue
continue

write(*,*) d2(1),resid(1),tang(1,1) -------------------------------------------!!!!!!! this write command recorded two times for one step

The values I got,
dis resid tang
1 100 100
1 100 100
2 200 100
2 200 100
3 300 100
3 300 100
.....

This has become a problem since when I was building my own element, I used state variables and at each step, I figured out it recorded two sets of resid and tang as below, only the second group is correct. However, both sets of values will be passed back to the main structure, leading to a failure of calculation. I wonder what could be the problems?

dis tang(1,1)
0 8438.91045339046
1 5530.69007988831
1 5530.69007988831
2 3712.52344293654 ! wrong
2 11595.7823440759 ! correct
3 3013.95024188580 ! wrong
3 11451.6296326108 ! correct
4 2440.38070368287 ! wrong
4 11409.7409248913 ! correct
....

[fmk]

are you using a linear solution algorithm?

[dray]

Hi Frank,

This is the code I used to build trussF.dll, I revised it to a 2D 3DOF element. I was using algorithm Newton, and I changed it to Linear now. But both gave me the results as I described.

=================
SUBROUTINE trussf(eleObj,modl,tang,resid,isw,error)

!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: TRUSSF
!DEC$ END IF

use elementTypes
use elementAPI
implicit none

type(eleObject)::eleObj
type(modelState)::modl
double precision tang(6,6)
double precision resid(6)
integer::isw;
integer::error;

integer :: tag, nd1, nd2, matTag, numCrd, i, j, numDOF
real *8, pointer::theParam(:)
integer, pointer::theNodes(:)

double precision A, dx, dy, L, cs, sn
double precision dLength, force, k

integer :: iData(3);
integer :: matTags(2);
c integer (c_int), target :: matTags(2);

type(c_ptr) :: theCMatPtr
type(c_ptr), pointer :: theCMatPtrPtr(:)
type(matObject), pointer :: theMat

double precision dData(1), nd1Crd(2), nd2Crd(2)
double precision d1(3), d2(3), tran(6)
double precision strs(1), strn(1), tng(1)


integer numData, err, matType

open(unit=111, file='truss.txt')

c outside functions called
c integer OPS_GetIntInput, OPS_GetDoubleInput, OPS_InvokeMaterial
c integer OPS_GetNodeCrd, OPS_AllocateElement, OPS_GetNodeDisp

IF (isw.eq.ISW_INIT) THEN

c get the input data - tag? nd1? nd2? A? matTag?

numData = 3
err = OPS_GetIntInput(numData, iData)
tag = iData(1);
nd1 = iData(2);
nd2 = iData(3);

numData = 1
err = OPS_GetDoubleInput(numData, dData)
A = dData(1);

numData = 1
err = OPS_GetIntInput(numData, iData)
matTag = iData(1);

c Allocate the element state

eleObj%tag = tag
eleObj%nnode = 2
eleObj%ndof = 6
eleObj%nparam = 4
eleObj%nstate = 0
eleObj%nmat = 1

matTags(1) = matTag;
matType = OPS_UNIAXIAL_MATERIAL_TYPE;
err = OPS_AllocateElement(eleObj, matTags, matType)

c theCMatPtr = theCMatPtrPtr(2);
c j=OPS_InvokeMaterialDirectly(theCMatPtr, modl, strn, strs,
c + tng, isw)

c element sets material functions
c call c_f_pointer(eleObj%mats, theCMatPtrPtr, [1]);
c theCMatPtrPtr(1) = theCMatPtr;

c Initialize the element properties

call c_f_pointer(eleObj%param, theParam, [4]);
call c_f_pointer(eleObj%node, theNodes, [2]);

numCrd = 2;
err = OPS_GetNodeCrd(nd1, numCrd, nd1Crd);
err = OPS_GetNodeCrd(nd2, numCrd, nd2Crd);

dx = nd2Crd(1)-nd1Crd(1);
dy = nd2Crd(2)-nd1Crd(2);

L = sqrt(dx*dx + dy*dy);

if (L == 0.0) then
c OPS_Error("Warning - truss element has zero length\n", 1);
return;
end if

cs = dx/L;
sn = dy/L;

theParam(1) = A;
theParam(2) = L;
theParam(3) = cs;
theParam(4) = sn;

theNodes(1) = nd1;
theNodes(2) = nd2;

ELSE

IF (isw == ISW_COMMIT) THEN

call c_f_pointer(eleObj%mats, theCMatPtrPtr, [1]);
theCMatPtr = theCMatPtrPtr(1);

j=OPS_InvokeMaterialDirectly(theCMatPtr, modl, strn, strs,
+ tng, isw)

ELSE IF (isw == ISW_REVERT_TO_START) THEN

call c_f_pointer(eleObj%mats, theCMatPtrPtr, [1]);
theCMatPtr = theCMatPtrPtr(1);

j=OPS_InvokeMaterialDirectly(theCMatPtr, modl, strn, strs,
+ tng, isw)

ELSE IF (isw == ISW_FORM_TANG_AND_RESID) THEN

call c_f_pointer(eleObj%param, theParam, [4]);
call c_f_pointer(eleObj%node, theNodes, [2]);
call c_f_pointer(eleObj%mats, theCMatPtrPtr, [1]);
theCMatPtr = theCMatPtrPtr(1);

A = theParam(1);
L = theParam(2);
cs = theParam(3);
sn = theParam(4);
nd1 = theNodes(1);
nd2 = theNodes(2);

numDOF = 3;
err = OPS_GetNodeDisp(nd1, numDOF, d1);
err = OPS_GetNodeDisp(nd2, numDOF, d2);

tran(1) = -cs;
tran(2) = -sn;
tran(3) = cs;
tran(4) = sn;

dLength = 0.0;
do 10 i = 1,2
dLength = dLength - (d2(i)-d1(i)) * tran(i);
10 continue

strn(1) = dLength/L;

c i = 0
c i=OPS_InvokeMaterial(eleObj, i, modl, strn, strs, tng, isw)
j=OPS_InvokeMaterialDirectly(theCMatPtr, modl, strn, strs,
+ tng, isw)

force = A*strs(1);
k = A*tng(1)/L;

resid(1)=force*tran(1);
resid(2)=force*tran(2);
resid(3)=0;
resid(4)=force*tran(3);
resid(5)=force*tran(4);
resid(6)=0;

tang(1,1)=k*tran(1)*tran(1);
tang(1,2)=k*tran(1)*tran(2);
tang(1,3)=0
tang(1,4)=k*tran(1)*tran(3);
tang(1,5)=k*tran(1)*tran(4);
tang(1,6)=0
tang(2,1)=k*tran(2)*tran(1);
tang(2,2)=k*tran(2)*tran(2);
tang(2,3)=0
tang(2,4)=k*tran(2)*tran(3);
tang(2,5)=k*tran(2)*tran(4);
tang(2,6)=0
tang(3,1)=0
tang(3,2)=0
tang(3,3)=0
tang(3,4)=0
tang(3,5)=0
tang(3,6)=0
tang(4,1)=k*tran(1)*tran(3);
tang(4,2)=k*tran(2)*tran(3);
tang(4,3)=0
tang(4,4)=k*tran(3)*tran(3);
tang(4,5)=k*tran(3)*tran(4);
tang(4,6)=0
tang(5,1)=k*tran(1)*tran(4);
tang(5,2)=k*tran(2)*tran(4);
tang(5,3)=0
tang(5,4)=k*tran(3)*tran(4);
tang(5,5)=k*tran(4)*tran(4);
tang(5,6)=0
tang(6,1)=0
tang(6,2)=0
tang(6,3)=0
tang(6,4)=0
tang(6,5)=0
tang(6,6)=0


write(111,*) d2(1),resid(1),tang(1,1)


END IF

END IF

c return error code
error = 0

END SUBROUTINE trussf

==================

This is the tcl code I used to run the model.

#create the ModelBuilder object
model basic -ndm 2 -ndf 3

# build the model

# add nodes - command: node nodeId xCrd yCrd
node 1 0.0 0.0
node 2 100 0

# add material - command: material <matType> matID <matArgs>
uniaxialMaterial Elastic 1 3000

# add truss elements - command: truss trussID node1 node2 A matID
element trussf 1 1 2 10.0 1


# set the boundary conditions - command: fix nodeID xResrnt? yRestrnt?
fix 1 1 1 1
fix 2 0 1 1

pattern Plain 1 Linear {
# apply the load - command: load nodeID xForce yForce
load 2 100 0 0
}

# build the components for the analysis object
system BandGeneral
constraints Plain
numberer RCM
test NormDispIncr 1.0e-7 10 0
algorithm Linear (Newton)

# create a Recorder object for the nodal displacements at node 4
recorder Node -file Example.out -load -node 2 -dof 1 2 disp
recorder Node -file RFrc.out -time -node 1 -dof 1 2 reaction

puts "analyze"
# perform the analysis
integrator DisplacementControl 2 1 1
# create the analysis object
analysis Static
analyze 6
integrator DisplacementControl 2 1 -1
analyze 12
integrator DisplacementControl 2 1 1
analyze 18
integrator DisplacementControl 2 1 -1
analyze 12

wipe


====================
The file "truss.txt", each step is recorded twice.

0.000000000000000E+000 0.000000000000000E+000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
6.00000000000000 -1800.00000000000 300.000000000000
6.00000000000000 -1800.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
0.000000000000000E+000 0.000000000000000E+000 300.000000000000
0.000000000000000E+000 0.000000000000000E+000 300.000000000000
-1.00000000000000 300.000000000000 300.000000000000
-1.00000000000000 300.000000000000 300.000000000000
-2.00000000000000 600.000000000000 300.000000000000
-2.00000000000000 600.000000000000 300.000000000000
-3.00000000000000 900.000000000000 300.000000000000
-3.00000000000000 900.000000000000 300.000000000000
-4.00000000000000 1200.00000000000 300.000000000000
-4.00000000000000 1200.00000000000 300.000000000000
-5.00000000000000 1500.00000000000 300.000000000000
-5.00000000000000 1500.00000000000 300.000000000000
-6.00000000000000 1800.00000000000 300.000000000000
-6.00000000000000 1800.00000000000 300.000000000000
-5.00000000000000 1500.00000000000 300.000000000000
-5.00000000000000 1500.00000000000 300.000000000000
-4.00000000000000 1200.00000000000 300.000000000000
-4.00000000000000 1200.00000000000 300.000000000000
-3.00000000000000 900.000000000000 300.000000000000
-3.00000000000000 900.000000000000 300.000000000000
-2.00000000000000 600.000000000000 300.000000000000
-2.00000000000000 600.000000000000 300.000000000000
-1.00000000000000 300.000000000000 300.000000000000
-1.00000000000000 300.000000000000 300.000000000000
0.000000000000000E+000 0.000000000000000E+000 300.000000000000
0.000000000000000E+000 0.000000000000000E+000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
6.00000000000000 -1800.00000000000 300.000000000000
6.00000000000000 -1800.00000000000 300.000000000000
7.00000000000000 -2100.00000000000 300.000000000000
7.00000000000000 -2100.00000000000 300.000000000000
8.00000000000000 -2400.00000000000 300.000000000000
8.00000000000000 -2400.00000000000 300.000000000000
9.00000000000000 -2700.00000000000 300.000000000000
9.00000000000000 -2700.00000000000 300.000000000000
10.0000000000000 -3000.00000000000 300.000000000000
10.0000000000000 -3000.00000000000 300.000000000000
11.0000000000000 -3300.00000000000 300.000000000000
11.0000000000000 -3300.00000000000 300.000000000000
12.0000000000000 -3600.00000000000 300.000000000000
12.0000000000000 -3600.00000000000 300.000000000000
11.0000000000000 -3300.00000000000 300.000000000000
11.0000000000000 -3300.00000000000 300.000000000000
10.0000000000000 -3000.00000000000 300.000000000000
10.0000000000000 -3000.00000000000 300.000000000000
9.00000000000000 -2700.00000000000 300.000000000000
9.00000000000000 -2700.00000000000 300.000000000000
8.00000000000000 -2400.00000000000 300.000000000000
8.00000000000000 -2400.00000000000 300.000000000000
7.00000000000000 -2100.00000000000 300.000000000000
7.00000000000000 -2100.00000000000 300.000000000000
6.00000000000000 -1800.00000000000 300.000000000000
6.00000000000000 -1800.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
5.00000000000000 -1500.00000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
4.00000000000000 -1200.00000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
3.00000000000000 -900.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
2.00000000000000 -600.000000000000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
1.00000000000000 -300.000000000000 300.000000000000
0.000000000000000E+000 0.000000000000000E+000 300.000000000000
0.000000000000000E+000 0.000000000000000E+000 300.000000000000

[fmk]

the problem is newStep() in the solution algoritth is causing the method to be invoked and then it is getting called in the solverCurrentStep of the algo .. a solution would be to store the current disp as a variable and then only do the code if the displ changes