Geometric Transformation of New Element in Fortran

[dray]

Hi all,

I`ve successfully compiled a two-node element in Fortran for OpenSees. The element aims to get the node displacements from OpenSees, calculate the stiffness matrix and residual forces, and then transfer the values back to OpenSees. The code is listed as below. When I run a simple analysis and apply the displacement control on one node. The results are perfect. However, when I put this element in a structure and apply displacement control on any other nodes. The results become unreasonable. What could be the possible reasons? I think one might be that there`s no geometric transformation between the natural coordinate and the global one. How can I get that part in the fortran code? I wonder if I can get some help here. Thank you so much.

In OpenSees, the element is called in the following code:
element abst 1 1 2 20 1
number node1 node2 area material

in which area and material can be ignored. All values are calulated in "UEL" function by given the node displacement of 1 and 2.

SUBROUTINE abst(eleObj,modl,tang,resid,isw,error)

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

use elementTypes
use elementAPI
implicit none

type(eleObject)::eleObj
type(modelState)::modl
double precision tang(4, 4)
double precision resid(4)
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(2), d2(2), tran(4)
double precision strs(1), strn(1), tng(1)

integer numData, err, matType
INTEGER ii,jj
DOUBLE PRECISION UU(4)
INTEGER iii

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 = 4
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
ELSE IF (isw == ISW_REVERT_TO_START) THEN
ELSE IF (isw == ISW_FORM_TANG_AND_RESID) THEN
call c_f_pointer(eleObj%node, theNodes, [2]);
nd1 = theNodes(1);
nd2 = theNodes(2);

call UEL(resid,tang,UU)

numDOF = 2;
err = OPS_GetNodeDisp(nd1, numDOF, d1);
err = OPS_GetNodeDisp(nd2, numDOF, d2);
UU(1)=d1(1)
UU(2)=d1(2)
UU(3)=d2(1)
UU(4)=d2(2)
END IF

END IF

c return error code
error = 0

END SUBROUTINE abst

SUBROUTINE UEL(RHS,AMATRX,U)
...

[fmk]

did you try a few of your elements and no other type? .. if that works, check the sign of the resisting force