int
FiberSection2d::addFiber(Fiber &newFiber)
{
  // need to create larger arrays
  int newSize = numFibers+1;
  UniaxialMaterial **newArray = new UniaxialMaterial *[newSize]; 
  double *newMatData = new double [2 * newSize];
  if (newArray == 0 || newMatData == 0) {
    opserr <<"FiberSection2d::addFiber -- failed to allocate Fiber pointers\n";
    return -1;
  }

  // copy the old pointers and data
  int i;
  for (i = 0; i < numFibers; i++) {
    newArray[i] = theMaterials[i];
    newMatData[2*i] = matData[2*i];
    newMatData[2*i+1] = matData[2*i+1];
  }

  // set the new pointers and data
  double yLoc, zLoc, Area;
  newFiber.getFiberLocation(yLoc, zLoc);
  Area = newFiber.getArea();
  newMatData[numFibers*2] = yLoc;
  newMatData[numFibers*2+1] = Area;
  UniaxialMaterial *theMat = newFiber.getMaterial();
  newArray[numFibers] = theMat->getCopy();

  if (newArray[numFibers] == 0) {
    opserr <<"FiberSection2d::addFiber -- failed to get copy of a Material\n";
    delete [] newMatData;
    return -1;
  }

  numFibers++;

  if (theMaterials != 0) {
    delete [] theMaterials;
    delete [] matData;
  }

  theMaterials = newArray;
  matData = newMatData;

  double Qz = 0.0;
  double A  = 0.0;

  // Recompute centroid
  for (i = 0; i < numFibers; i++) {
    yLoc = -matData[2*i];
    Area = matData[2*i+1];
    A  += Area;
    Qz += yLoc*Area;
  }

  yBar = Qz/A;

  return 0;
}

// constructor:
//  responsible for allocating the necessary space needed by each object
//  and storing the tags of the CorotTruss end nodes.
CorotTruss::CorotTruss(int tag, int dim,
			   int Nd1, int Nd2, 
			   UniaxialMaterial &theMat,
			   double a, double r)
  :Element(tag,ELE_TAG_CorotTruss),     
  theMaterial(0), connectedExternalNodes(2),
  numDOF(0), numDIM(dim),
  Lo(0.0), Ln(0.0), 
  A(a), rho(r), R(3,3),
  theMatrix(0), theVector(0)
{
  // get a copy of the material and check we obtained a valid copy
  theMaterial = theMat.getCopy();
  if (theMaterial == 0) {
    opserr << "FATAL CorotTruss::CorotTruss - " <<  tag <<
      "failed to get a copy of material with tag " << theMat.getTag() << endln;
    exit(-1);
  }
  
  // ensure the connectedExternalNode ID is of correct size & set values
  if (connectedExternalNodes.Size() != 2) {
    opserr << "FATAL CorotTruss::CorotTruss - " <<  tag <<
      "failed to create an ID of size 2\n";
    exit(-1);
  }
  
  connectedExternalNodes(0) = Nd1;
  connectedExternalNodes(1) = Nd2;        

  // set node pointers to NULL
  theNodes[0] = 0;
  theNodes[1] = 0;
}

matObj *OPS_GetMaterial(int *matTag, int *matType)
{

  if (*matType == OPS_UNIAXIAL_MATERIAL_TYPE) {
    UniaxialMaterial *theUniaxialMaterial = OPS_getUniaxialMaterial(*matTag);
    
    if (theUniaxialMaterial != 0) {
      
      UniaxialMaterial *theCopy = theUniaxialMaterial->getCopy();
      //  uniaxialMaterialObjectCount++;
      // theUniaxialMaterials[uniaxialMaterialObjectCount] = theCopy;
      
      matObject *theMatObject = new matObject;
      theMatObject->tag = *matTag;
      theMatObject->nParam = 1;
      theMatObject->nState = 0;
      
      theMatObject->theParam = new double[1];
      //  theMatObject->theParam[0] = uniaxialMaterialObjectCount;
      theMatObject->theParam[0] = 1; // code for uniaxial material
      
      theMatObject->tState = 0;
      theMatObject->cState = 0;
      theMatObject->matFunctPtr = OPS_InvokeMaterialObject;

      theMatObject->matObjectPtr = theCopy;
      
      return theMatObject;
    }
    
    fprintf(stderr,"getMaterial - no uniaxial material exists with tag %d\n", *matTag);    
    return 0;

  } else if (*matType == OPS_SECTION_TYPE) {
    fprintf(stderr,"getMaterial - not yet implemented for Section\n");    
    return 0;
  } else {

    //    NDMaterial *theNDMaterial = theModelBuilder->getNDMaterial(*matTag);

    //    if (theNDMaterial != 0) 
      //      theNDMaterial = theNDMaterial->getCopy(matType);
      //    else {
      //      fprintf(stderr,"getMaterial - no nd material exists with tag %d\n", *matTag);          
      //      return 0;
      //    }

      //    if (theNDMaterial == 0) {
    //      fprintf(stderr,"getMaterial - material with tag %d cannot deal with %d\n", *matTag, matType);          
    //      return 0;
    //    }

    fprintf(stderr,"getMaterial - not yet implemented for nDMaterial\n");    
    return 0;
  }

  fprintf(stderr,"getMaterial - unknown material type\n");    
  return 0;

}

// constructors:
FiberSection3d::FiberSection3d(int tag, int num, Fiber **fibers): 
  SectionForceDeformation(tag, SEC_TAG_FiberSection3d),
  numFibers(num), sizeFibers(num), theMaterials(0), matData(0),
  QzBar(0.0), QyBar(0.0), Abar(0.0), yBar(0.0), zBar(0.0), sectionIntegr(0), e(3), s(0), ks(0)
{
  if (numFibers != 0) {
    theMaterials = new UniaxialMaterial *[numFibers];

    if (theMaterials == 0) {
      opserr << "FiberSection3d::FiberSection3d -- failed to allocate Material pointers\n";
      exit(-1);
    }

    matData = new double [numFibers*3];

    if (matData == 0) {
      opserr << "FiberSection3d::FiberSection3d -- failed to allocate double array for material data\n";
      exit(-1);
    }

    for (int i = 0; i < numFibers; i++) {
      Fiber *theFiber = fibers[i];
      double yLoc, zLoc, Area;
      theFiber->getFiberLocation(yLoc, zLoc);
      Area = theFiber->getArea();

      QzBar += yLoc*Area;
      QyBar += zLoc*Area;
      Abar  += Area;

      matData[i*3] = yLoc;
      matData[i*3+1] = zLoc;
      matData[i*3+2] = Area;
      UniaxialMaterial *theMat = theFiber->getMaterial();
      theMaterials[i] = theMat->getCopy();

      if (theMaterials[i] == 0) {
	opserr << "FiberSection3d::FiberSection3d -- failed to get copy of a Material\n";
	exit(-1);
      }
    }

    yBar = QzBar/Abar;
    zBar = QyBar/Abar;
  }

  s = new Vector(sData, 3);
  ks = new Matrix(kData, 3, 3);

  sData[0] = 0.0;
  sData[1] = 0.0;
  sData[2] = 0.0;

  for (int i=0; i<9; i++)
    kData[i] = 0.0;

  code(0) = SECTION_RESPONSE_P;
  code(1) = SECTION_RESPONSE_MZ;
  code(2) = SECTION_RESPONSE_MY;
}

static 
void OPS_InvokeMaterialObject(struct matObject *theMat, modelState *theModel,double *strain, double *tang, double *stress, int *isw, int *result)
{
  int matType = theMat->theParam[0];

  if (matType == 1) {
    //  UniaxialMaterial *theMaterial = theUniaxialMaterials[matCount];
    UniaxialMaterial *theMaterial = (UniaxialMaterial *)theMat->matObjectPtr;
    if (theMaterial == 0) {
      *result = -1;
      return;
    }
    
    if (*isw == ISW_COMMIT) {
      *result =  theMaterial->commitState();
      return;
    } else if (*isw == ISW_REVERT) {
      *result =  theMaterial->revertToLastCommit();
      return;
    } else if (*isw == ISW_REVERT_TO_START) {
      *result =  theMaterial->revertToStart();
      return;
    } else if (*isw == ISW_FORM_TANG_AND_RESID) {
      double matStress = 0.0;
      double matTangent = 0.0;
      int res = theMaterial->setTrial(strain[0], matStress, matTangent);
      stress[0] = matStress;
      tang[0] = matTangent;
      *result = res;
      return;
    }
  }
  
  return;
}

int OPS_setStrain()
{
    if (OPS_GetNumRemainingInputArgs() != 1) {
	opserr<<"testUniaxialMaterial - You must provide a strain value.\n";
	return -1;
    }

    UniaxialMaterial* material = theTestingUniaxialMaterial;

    if (material == 0) {
	opserr<<"setStrain WARNING no active UniaxialMaterial - use testUniaxialMaterial command.\n";
	return -1;
    }

    double strain;
    int numData = 1;
    if (OPS_GetDoubleInput(&numData, &strain) < 0) {
	opserr<<"invalid double value\n";
	return -1;
    }

    material->setTrialStrain(strain);
    material->commitState();

    return 0;
}

int
FiberSection3dThermal::commitSensitivity(const Vector& defSens, int gradIndex, int numGrads)
{

  // here add SHVs to store the strain sensitivity.

  if (SHVs == 0) {
    SHVs = new Matrix(3,numGrads);
  }
  
  (*SHVs)(0,gradIndex) = defSens(0);
  (*SHVs)(1,gradIndex) = defSens(1);
  (*SHVs)(2,gradIndex) = defSens(2);

  int loc = 0;

  double d0 = defSens(0);
  double d1 = defSens(1);
  double d2 = defSens(2);

  for (int i = 0; i < numFibers; i++) {
    UniaxialMaterial *theMat = theMaterials[i];
   	double y = matData[loc++] - yBar;
	double z = matData[loc++] - zBar;
	loc++;   // skip A data.

	double strainSens = d0 + y*d1 + z*d2;


    
	theMat->commitSensitivity(strainSens,gradIndex,numGrads);
  }

  return 0;
}

// constructors:
FiberSection2d::FiberSection2d(int tag, int num, Fiber **fibers): 
  SectionForceDeformation(tag, SEC_TAG_FiberSection2d),
  numFibers(num), theMaterials(0), matData(0),
  yBar(0.0), sectionIntegr(0), e(2), eCommit(2), s(0), ks(0), dedh(2)
{
  if (numFibers != 0) {
    theMaterials = new UniaxialMaterial *[numFibers];

    if (theMaterials == 0) {
      opserr << "FiberSection2d::FiberSection2d -- failed to allocate Material pointers";
      exit(-1);
    }

    matData = new double [numFibers*2];

    if (matData == 0) {
      opserr << "FiberSection2d::FiberSection2d -- failed to allocate double array for material data\n";
      exit(-1);
    }


    double Qz = 0.0;
    double A  = 0.0;
    
    for (int i = 0; i < numFibers; i++) {
      Fiber *theFiber = fibers[i];
      double yLoc, zLoc, Area;
      theFiber->getFiberLocation(yLoc, zLoc);
      Area = theFiber->getArea();
      A  += Area;
      Qz += yLoc*Area;
      matData[i*2] = yLoc;
      matData[i*2+1] = Area;
      UniaxialMaterial *theMat = theFiber->getMaterial();
      theMaterials[i] = theMat->getCopy();

      if (theMaterials[i] == 0) {
	opserr << "FiberSection2d::FiberSection2d -- failed to get copy of a Material\n";
	exit(-1);
      }
    }    

    yBar = Qz/A;  
  }

  s = new Vector(sData, 2);
  ks = new Matrix(kData, 2, 2);

  sData[0] = 0.0;
  sData[1] = 0.0;

  kData[0] = 0.0;
  kData[1] = 0.0;
  kData[2] = 0.0;
  kData[3] = 0.0;

  code(0) = SECTION_RESPONSE_P;
  code(1) = SECTION_RESPONSE_MZ;
}

const Matrix&
FiberSection3d::getInitialTangent(void)
{
  static double kInitialData[9];
  static Matrix kInitial(kInitialData, 3, 3);
  
  kInitialData[0] = 0.0; kInitialData[1] = 0.0; 
  kInitialData[2] = 0.0; kInitialData[3] = 0.0;
  kInitialData[4] = 0.0; kInitialData[5] = 0.0; 
  kInitialData[6] = 0.0; kInitialData[7] = 0.0;
  kInitialData[8] = 0.0; 

  static double yLocs[10000];
  static double zLocs[10000];
  static double fiberArea[10000];

  if (sectionIntegr != 0) {
    sectionIntegr->getFiberLocations(numFibers, yLocs, zLocs);
    sectionIntegr->getFiberWeights(numFibers, fiberArea);
  }  
  else {
    for (int i = 0; i < numFibers; i++) {
      yLocs[i] = matData[3*i];
      zLocs[i] = matData[3*i+1];
      fiberArea[i] = matData[3*i+2];
    }
  }

  for (int i = 0; i < numFibers; i++) {
    UniaxialMaterial *theMat = theMaterials[i];
    double y = yLocs[i] - yBar;
    double z = zLocs[i] - zBar;
    double A = fiberArea[i];

    double tangent = theMat->getInitialTangent();

    double value = tangent * A;
    double vas1 = -y*value;
    double vas2 = z*value;
    double vas1as2 = vas1*z;

    kInitialData[0] += value;
    kInitialData[1] += vas1;
    kInitialData[2] += vas2;
    
    kInitialData[4] += vas1 * -y;
    kInitialData[5] += vas1as2;
    
    kInitialData[8] += vas2 * z; 
  }

  kInitialData[3] = kInitialData[1];
  kInitialData[6] = kInitialData[2];
  kInitialData[7] = kInitialData[5];

  return kInitial;
}

int
FiberSection2d::sendSelf(int commitTag, Channel &theChannel)
{
  int res = 0;

  // create an id to send objects tag and numFibers, 
  //     size 3 so no conflict with matData below if just 1 fiber
  static ID data(3);
  data(0) = this->getTag();
  data(1) = numFibers;
  int dbTag = this->getDbTag();
  res += theChannel.sendID(dbTag, commitTag, data);
  if (res < 0) {
    opserr <<  "FiberSection2d::sendSelf - failed to send ID data\n";
    return res;
  }    

  if (numFibers != 0) {
    
    // create an id containingg classTag and dbTag for each material & send it
    ID materialData(2*numFibers);
    for (int i=0; i<numFibers; i++) {
      UniaxialMaterial *theMat = theMaterials[i];
      materialData(2*i) = theMat->getClassTag();
      int matDbTag = theMat->getDbTag();
      if (matDbTag == 0) {
	matDbTag = theChannel.getDbTag();
	if (matDbTag != 0)
	  theMat->setDbTag(matDbTag);
      }
      materialData(2*i+1) = matDbTag;
    }    
    
    res += theChannel.sendID(dbTag, commitTag, materialData);
    if (res < 0) {
      opserr <<  "FiberSection2d::sendSelf - failed to send material data\n";
      return res;
    }    

    // send the fiber data, i.e. area and loc
    Vector fiberData(matData, 2*numFibers);
    res += theChannel.sendVector(dbTag, commitTag, fiberData);
    if (res < 0) {
      opserr <<  "FiberSection2d::sendSelf - failed to send material data\n";
      return res;
    }    

    // now invoke send(0 on all the materials
    for (int j=0; j<numFibers; j++)
      theMaterials[j]->sendSelf(commitTag, theChannel);

  }

  return res;
}

int
FiberSection3dThermal::revertToLastCommit(void)
{
  int err = 0;

  // Last committed section deformations
  e = eCommit;


  kData[0] = 0.0; kData[1] = 0.0; kData[2] = 0.0; kData[3] = 0.0;
  kData[4] = 0.0; kData[5] = 0.0; kData[6] = 0.0; kData[7] = 0.0;
  kData[8] = 0.0; 
  sData[0] = 0.0; sData[1] = 0.0;  sData[2] = 0.0; 

  int loc = 0;

  for (int i = 0; i < numFibers; i++) {
    UniaxialMaterial *theMat = theMaterials[i];
    double y = matData[loc++] - yBar;
    double z = matData[loc++] - zBar;
    double A = matData[loc++];

    // invoke revertToLast on the material
    err += theMat->revertToLastCommit();

    double tangent = theMat->getTangent();
    double stress = theMat->getStress();

    double value = tangent * A;
    double vas1 = y*value;
    double vas2 = z*value;
    double vas1as2 = vas1*z;

    kData[0] += value;
    kData[1] += vas1;
    kData[2] += vas2;
    
    kData[4] += vas1 * y;
    kData[5] += vas1as2;
    
    kData[8] += vas2 * z; 

    double fs0 = stress * A;
    sData[0] += fs0;
    sData[1] += fs0 * y;
    sData[2] += fs0 * z;
  }

  kData[3] = kData[1];
  kData[6] = kData[2];
  kData[7] = kData[5];

  return err;
}

int
FiberSection2d::setTrialSectionDeformation (const Vector &deforms)
{
  int res = 0;

  e = deforms;

  kData[0] = 0.0; kData[1] = 0.0; kData[2] = 0.0; kData[3] = 0.0;
  sData[0] = 0.0; sData[1] = 0.0;

  double d0 = deforms(0);
  double d1 = deforms(1);

  static double fiberLocs[10000];
  static double fiberArea[10000];

  if (sectionIntegr != 0) {
    sectionIntegr->getFiberLocations(numFibers, fiberLocs);
    sectionIntegr->getFiberWeights(numFibers, fiberArea);
  }  
  else {
    for (int i = 0; i < numFibers; i++) {
      fiberLocs[i] = matData[2*i];
      fiberArea[i] = matData[2*i+1];
    }
  }
  
  for (int i = 0; i < numFibers; i++) {
    UniaxialMaterial *theMat = theMaterials[i];
    double y = fiberLocs[i] - yBar;
    double A = fiberArea[i];

    // determine material strain and set it
    double strain = d0 - y*d1;
    double tangent, stress;
    res += theMat->setTrial(strain, stress, tangent);

    double ks0 = tangent * A;
    double ks1 = ks0 * -y;
    kData[0] += ks0;
    kData[1] += ks1;
    kData[3] += ks1 * -y;

    double fs0 = stress * A;
    sData[0] += fs0;
    sData[1] += fs0 * -y;
  }

  kData[2] = kData[1];

  return res;
}

int
FiberSection3d::commitSensitivity(const Vector& defSens, int gradIndex, int numGrads)
{

  double d0 = defSens(0);
  double d1 = defSens(1);
  double d2 = defSens(2);

  //dedh = defSens;

  static double yLocs[10000];
  static double zLocs[10000];

  if (sectionIntegr != 0)
    sectionIntegr->getFiberLocations(numFibers, yLocs, zLocs);
  else {
    for (int i = 0; i < numFibers; i++) {
      yLocs[i] = matData[3*i];
      zLocs[i] = matData[3*i+1];
    }
  }

  static double dydh[10000];
  static double dzdh[10000];

  if (sectionIntegr != 0)
    sectionIntegr->getLocationsDeriv(numFibers, dydh, dzdh);  
  else {
    for (int i = 0; i < numFibers; i++) {
      dydh[i] = 0.0;
      dzdh[i] = 0.0;
    }
  }

  double y, z;

  double depsdh = 0;

  for (int i = 0; i < numFibers; i++) {
    UniaxialMaterial *theMat = theMaterials[i];
    y = yLocs[i] - yBar;
    z = zLocs[i] - zBar;

    // determine material strain and set it
    depsdh = d0 - y*d1 + z*d2 - dydh[i]*e(1) + dzdh[i]*e(2);

    theMat->commitSensitivity(depsdh,gradIndex,numGrads);
  }

  return 0;
}

// constructor:
//  responsible for allocating the necessary space needed by each object
//  and storing the tags of the CorotTruss end nodes.
XC::CorotTruss::CorotTruss(int tag, int dim,int Nd1, int Nd2, UniaxialMaterial &theMat,double a)
  :CorotTrussBase(tag,ELE_TAG_CorotTruss,dim,Nd1,Nd2), theMaterial(nullptr), A(a)
  {
    // get a copy of the material and check we obtained a valid copy
    theMaterial = theMat.getCopy();
    if(theMaterial == 0)
      {
        std::cerr << getClassName() << "::" << __FUNCTION__
		  << "; FATAL error in element: " <<  tag
		  << "failed to get a copy of material with tag "
		  << theMat.getTag() << std::endl;
        exit(-1);
      }
  }

Truss2::Truss2(int tag, 
	int dim,
	int Nd1, int Nd2, int oNd1, int oNd2, 
	UniaxialMaterial &theMat,
	double a, double r, int damp)
	:Element(tag,ELE_TAG_Truss2),     
	theMaterial(0), theBetaMaterial(0), connectedExternalNodes(2), connectedExternalOtherNodes(2),
	dimension(dim), numDOF(0), theLoad(0),
	theMatrix(0), theVector(0),
	L(0.0), A(a), rho(r), doRayleighDamping(damp)
{
	// get a copy of the material and check we obtained a valid copy
	theMaterial = theMat.getCopy();
	if (theMaterial == 0) {
	  opserr << "FATAL Truss2::Truss2 - " << tag <<
	    "failed to get a copy of material with tag " << theMat.getTag() << endln;
	  exit(-1);
	} else if (theMaterial->getClassTag() == MAT_TAG_ConcretewBeta) {
	  theBetaMaterial = (ConcretewBeta *) theMaterial;
	}

	// ensure the connectedExternalNode ID is of correct size & set values
	if (connectedExternalNodes.Size() != 2 || connectedExternalOtherNodes.Size() != 2) {
		opserr << "FATAL Truss2::Truss2 - " <<  tag << "failed to create an ID of size 2\n";
		exit(-1);
	}

	connectedExternalNodes(0) = Nd1;
	connectedExternalNodes(1) = Nd2; 

	/// some work to save the other nodes:
	connectedExternalOtherNodes(0) = oNd1;
	connectedExternalOtherNodes(1) = oNd2;

	// set node pointers to NULL
	for (int i=0; i<2; i++) {
		theNodes[i] = 0;
		theOtherNodes[i] = 0;
	}

	cosX[0] = 0.0;
	cosX[1] = 0.0;
	cosX[2] = 0.0;

	// AddingSensitivity:BEGIN /////////////////////////////////////
	parameterID = 0;
	theLoadSens = 0;
	// AddingSensitivity:END //////////////////////////////////////
}

//  Construct element with one unidirectional material (numMaterials1d=1)
CoupledZeroLength::CoupledZeroLength(int tag,
			 int Nd1, int Nd2, 
			 UniaxialMaterial &theMat,
			 int direction1, int direction2,
			 int doRayleigh)
  :Element(tag,ELE_TAG_CoupledZeroLength),     
   connectedExternalNodes(2),
   dimension(0), numDOF(0), transformation(3,3), useRayleighDamping(doRayleigh),
   theMatrix(0), theVector(0),
   theMaterial(0), dirn1(direction1), dirn2(direction2), d0(0), v0(0)
{
  // allocate memory for numMaterials1d uniaxial material models
  theMaterial = theMat.getCopy();
  if ( theMaterial == 0) {
    opserr << "FATAL CoupledZeroLength::CoupledZeroLength - failed to create a 1d  material\n";
    exit(-1);
  }

  // initialize uniaxial materials and directions and check for valid values
  if (direction1 < 0 || direction1 > 5 || direction2 < 0 || direction2 > 5) {
    opserr << "FATAL: CoupledZeroLength::CoupledZeroLength - invalid diection\n";
    exit(-1);
  }

  connectedExternalNodes(0) = Nd1;
  connectedExternalNodes(1) = Nd2;
  dX = 0.0;
  dY = 0.0;
  fX = 0.0;
  fY = 0.0;
}

int
FiberSection2d::commitSensitivity(const Vector& defSens,
				  int gradIndex, int numGrads)
{
  double d0 = defSens(0);
  double d1 = defSens(1);

  dedh = defSens;

  static double fiberLocs[10000];

  if (sectionIntegr != 0)
    sectionIntegr->getFiberLocations(numFibers, fiberLocs);
  else {
    for (int i = 0; i < numFibers; i++)
      fiberLocs[i] = matData[2*i];
  }

  static double locsDeriv[10000];
  static double areaDeriv[10000];

  if (sectionIntegr != 0) {
    sectionIntegr->getLocationsDeriv(numFibers, locsDeriv);  
    sectionIntegr->getWeightsDeriv(numFibers, areaDeriv);
  }
  else {
    for (int i = 0; i < numFibers; i++) {
      locsDeriv[i] = 0.0;
      areaDeriv[i] = 0.0;
    }
  }

  double y;
  double kappa = e(1);

  for (int i = 0; i < numFibers; i++) {
    UniaxialMaterial *theMat = theMaterials[i];
    y = fiberLocs[i] - yBar;

    // determine material strain and set it
    double strainSens = d0 - y*d1 - locsDeriv[i]*kappa;
    theMat->commitSensitivity(strainSens,gradIndex,numGrads);
  }

  return 0;
}

N4BiaxialTruss::N4BiaxialTruss(int tag, 
int dim,
int Nd1, int Nd2, 
int GNd1, int GNd2, 
UniaxialMaterial &theMat,
double a, double r, int damp)
:Element(tag,ELE_TAG_N4BiaxialTruss),     
theMaterial_1(0), theBetaMaterial_1(0), 
theMaterial_2(0), theBetaMaterial_2(0), 
connectedExternalNodes(4),
dimension(dim), numDOF(0), theLoad(0),
theMatrix(0), theVector(0), theVector2(0),
L(0.0), A(a), rho(r), doRayleighDamping(damp)
{
	// get a copy of the material and check we obtained a valid copy
	theMaterial_1 = theMat.getCopy();
	theMaterial_2 = theMat.getCopy();
	if ((theMaterial_1 == 0) || (theMaterial_2 == 0)) {
		opserr << "FATAL N4BiaxialTruss::N4BiaxialTruss - " << tag <<
		"failed to get a copy of material with tag " << theMat.getTag() << endln;
		exit(-1);
	} else if (theMat.getClassTag() == MAT_TAG_ConcretewBeta) {
	  theBetaMaterial_1 = (ConcretewBeta *) theMaterial_1;
	  theBetaMaterial_2 = (ConcretewBeta *) theMaterial_2;
	}
	
	// ensure the connectedExternalNode ID is of correct size & set values
	if (connectedExternalNodes.Size() != 4) {
		opserr << "FATAL N4BiaxialTruss::N4BiaxialTruss - " <<  tag << "failed to create an node ID array of size 4\n";
		exit(-1);
	}

	connectedExternalNodes(0) = Nd1;
	connectedExternalNodes(1) = Nd2;
	connectedExternalNodes(2) = GNd1;
	connectedExternalNodes(3) = GNd2;

	// set node pointers to NULL
	for (int i=0; i<4; i++)
	theNodes[i] = 0;

	cosX[0] = 0.0;
	cosX[1] = 0.0;
	cosX[2] = 0.0;
}

void OPS_printUniaxialMaterial(OPS_Stream &s, int flag) {
  if (flag == OPS_PRINT_PRINTMODEL_JSON) {
    s << "\t\t\"uniaxialMaterials\": [\n";        
    MapOfTaggedObjectsIter theObjects = theUniaxialMaterialObjects.getIter();
    theObjects.reset();
    TaggedObject *theObject;
    int count = 0;
    int numComponents = theUniaxialMaterialObjects.getNumComponents();    
    while ((theObject = theObjects()) != 0) {
      UniaxialMaterial *theMaterial = (UniaxialMaterial *)theObject;
      theMaterial->Print(s, flag);
      if (count < numComponents-1)
	s << ",\n";
      count++;      
    }
    s << "\n\t\t]";
  }
}

int OPS_getStress()
{
    UniaxialMaterial* material = theTestingUniaxialMaterial;
    if (material == 0) {
	opserr<<"getStrain WARNING no active UniaxialMaterial - use testUniaxialMaterial command.\n";
	return -1;
    }

    double stress = material->getStress();

    int numData = 1;

    if (OPS_SetDoubleOutput(&numData, &stress) < 0) {
	opserr<<"failed to set stress\n";
	return -1;
    }

    return 0;
}