Adding Your Code to OpenSees
Frank McKenna
UC Berkeley
OpenSeesDays 2012
Outline
Introduction
Adding a New Material to [Link]
Adding a New Integrator to [Link]
Summary & Conclusions
NOTE: I HOPE NOT TO GET BOGGED
DOWN IN C++ ISSUES.
Traditional Finite Element
Software
Commercial fe codes are large complex
software containing over one million lines of code.
They are closed-source but do allow new
element and material routines to be added. (at
least the serious ones do)
They are slow to change as hardware changes
and they do not allow researchers to play with
other equally important aspects of the code.
They do not promote the active sharing of new
code.
THE ADVANTAGE OF A
SOFTWARE FRAMEWORK
SUCH AS OPENSEES IS THAT
YOU DONT HAVE TO
UNDERSTAND ALL OF IT TO
BUILD APPLICATIONS OR
MAKE CONTRIBUTIONS
YOU JUST NEED TO
UNDERSTAND THAT PART
THAT CONCERNS YOU
OPENSEES FEM ABSTRACT
CLASSES
Framework Contains both Abstract
Classes & Concrete Subclasses
Material
Uniaxial
Elastic
ElasticPP
Hardening
Concrete
Steel
Hysteretic
PY-TZ-QZ
Parallel
Series
Gap
Fatigue
nD
Element
Section
Elastic
Fiber
Elastic
J2
DruckerPrager
FluidSolidPorous
PressureMultiYield(dependent, independent)
(over 250 material classes)
Truss
ZeroLength
ElasticBeamColumn
ForceBeamColumn
DispBeamColumn
BeamWithHinges
Quad (many forms)
Shell
Brick(many forms)
Joint
Contact
Integrator
StaticIntegrator
LoadControl
DispControl
ArcLength
TransientIntegrator
CentralDifference
Newmark
HHT
GeneralizedAlhpa
NewmarkExplicit
TRBDF2
AlphaOS
(>100 element classes)
(>35 classes)
Unknown Class Type
When [Link] is running and it
comes across a class type it knows nothing
about BEFORE GIVING AN ERROR IT
WILL TRY AND LOAD A LIBRARY OF
THAT CLASSES NAME. If it cannot find
a library of the appropriate name or the
procedure of the appropriate name in the
library it will FAIL. Commands/Script
Output
To Add a New Class you must:
1) Provide Code that meets the Interface
of the appropriate super-class
2) you must BUILD THE LIBRARY
3) make it accessible to the program.
WHILE C++ IS THE GLUE LANGUAGE
THAT HOLDS OPENSEES TOGETHER
YOUR CODE DOES NOT HAVE TO BE
WRITTEN IN C++.
C and FORTRAN OPTIONS ARE ALSO
AVAILABLE
class UniaxialMaterial : public Material
{
public:
UniaxialMaterial(int tag, int classTag);
virtual ~UniaxialMaterial();
virtual
virtual
virtual
virtual
virtual
Must be overridden by
subclass, pure virtual
int setTrialStrain(double strain, double strainRate = 0.0) = 0;
double getStrain(void) = 0;
double getStress(void) = 0;
double getTangent(void) = 0;
double getInitialtangent(void) =0;
virtual int commitState(void) = 0;
virtual int revertToLastCommit(void) = 0;
virtual int revertToStart(void) = 0;
virtual UniaxialMaterial *getCopy(void) = 0;
Can be overridden by subclass
virtual Response *setResponse(const char **argv, int argc,OPS_Stream &theOutput);
virtual int getResponse(int responseID, Information &info);
virtual void Print(OPS_Stream &s, int flag =0);
virtual int sendSelf(int commitTag, Channel &theChannel)=0;
virtual int recvSelf(int commitTag, Chanel &theChannel, FEM_ObjectBroker &theBroker)=0;
// THERE ARE SOME OTHERS .. BUT THESE ARE PURE VIRTUAL ONES THAT MUST BE PROVIDED
protected:
private:
};
Adding New Code
For those new to Programming NEVER EVER NEVER START
WITH AN EMPTY FILE .. TAKE SOMETHING SIMILAR THAT
WORKS AND MAKE CHANGES TO THAT FILE.
We provide C++, C and Fortran examples on-line using svn
svn://[Link]/usr/local/svn/OpenSees/trunk/DEVELOPER
TortoiseSVN for Windows Users
your
Source Code Tree in
DEVELOPER
Makefile
WindowsInstructions
UnixInstructions
core
integrator
Trapezoidal.h
[Link]
recorder
cpp
material
c
element
fortran
ElasticPPcpp.h
elasticPPc.c
[Link] [Link]
[Link]
element
cpp
fortran
elasticPPf.f
[Link]
class ElasticPPcpp: public UniaxialMaterial {
public:
Material ElasticPPcpp(int tag, double e, double eyp);
ElasticPPcpp();
material input properties
Name
~ ElasticPPcpp();
int setTrialStrain(double strain, double strainRate=0.0);
double getStrain(void);
double getStress(void);
double getTangent(void);
double getInitialTangent(void);
int commitState(void);
int revertToLastCommit(void);
int revertToStart(void);
UniaxialMaterial *getCopy(void);
int sendSelf(int commitTag, Channel &theChannel);
int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker);
void Print(OPS_Stream &s, int flag = 0);
private:
data unique to
double fyp, fyn; // pos & neg yield stress
material includes:
double ezero, E,ep; // init strain, elastic mod
double ep; // plastic strain at last commit Material parameters,
double trialStrain, trialStress, trialTangent;
& State variables
double commitStrain, commitStress, commitTangent;
};
ElasticPPcpp::ElasticPPcpp(int tag, double e, double eyp)
:UniaxialMaterial(tag, 0),
ezero(0), E(e), ep(0.0) trialStrain(0.0),trialStress(0.0),trialTangent(e),
commitStreain(0.0),commitStress(0.0),commitTangent(e)
{
fyp=E*eyp;
fyn = -fyp;
}
ElasticPPcpp::ElasticPPcpp()
:UniaxialMaterial(tag, 0)
fyp(0),fyn(0),ezero(0), E(0),ep(0),
trialStrain(0.0),trialStress(0.0),trialTangent(0),
commitStrain(0.0),commitStress(0.0),commitTangent(e){
}
ElasticPPcpp::~ElasticPPcpp
{
// does nothing .. No memory to clean up
}
UniaxialMaterial *ElasticPPcpp::getCopy(void)
{
ElasticPPcpp *theCopy = new ElasticPPcpp(this->getTag(), E, fyp/E);
return theCopy;
};
Hardest Method to Write
ElasticPPcpp::setTrialStrain(double strain, double strainRate)
{
if (fabs(trialStrain - strain) < DBL_EPSILON)
return 0;
trialStrain = strain;
double sigtrial; // trial stress
double f;
// yield function
// compute trial stress
sigtrial = E * ( trialStrain - ezero - ep );
// evaluate yield function
if ( sigtrial >= 0.0 )
f = sigtrial - fyp;
else
f = -sigtrial + fyn;
double fYieldSurface = - E * DBL_EPSILON;
if ( f <= fYieldSurface ) {
// elastic
trialStress = sigtrial;
trialTangent = E;
} else {
// plastic
if ( sigtrial > 0.0 ) {
trialStress = fyp;
} else {
trialStress = fyn;
}
}
trialTangent = 0.0;
return 0;
double
ElasticPPcpp::getStrain(void)
{
return trialStrain;
}
double
ElasticPPcpp::getStress(void)
{
return trialStress;
}
double
ElasticPPcpp::getTangent(void)
{
return trialTangent;
}
int
ElasticPPcpp::revertToLastCommit(void)
{
trialStrain = commitStrain;
trialTangent = commitTangent;
trialStress = commitStress;
return 0;
}
Interpreter looking for this function in lib
OPS_Export void *
OPS_ElasticPPcpp()
{
You can give yourself some
if (numElasticPPcpp == 0) {
[Link]
reference
opserr << "ElasticPPcpp unaxialKUDOS,
material - Written
fmk UC Berkeley\n;
numElasticPPcpp =1;
if you used this
}
// Pointer to a uniaxial material that will be returned
UniaxialMaterial *theMaterial = 0;
int iData[1];
double dData[2];
parse the script for
int numData;
three material parameters
numData = 1;
if (OPS_GetIntInput(&numData, iData) != 0) {
opserr << "WARNING invalid uniaxialMaterial ElasticPP tag" << endln;
return 0;
}
numData = 2;
if (OPS_GetDoubleInput(&numData, dData) != 0) {
opserr << "WARNING invalid E & ep\n";
return 0;
}
Function returns new material
theMaterial = new ElasticPPcpp(iData[0], dData[0], dData[1]);
return theMaterial;
}
C & Fortran Procedural
Languages Can Also Be Used
OPS_Export void
elasticPPc (matObj *thisObj,
modelState *model,
double *strain,
double *tang,
double *stress,
int *isw,
int *result)
{
*result = 0;
if (*isw == ISW_INIT) {
SUBROUTINE ELASTICPPF(matObj,model,strain,tang,stress,isw,err
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: ELASTICPPF
!DEC$ END IF
use materialTypes
use materialAPI
implicit none
IF ([Link].ISW_INIT) THEN
c
get the input data - tag? E? eyp?
double dData[2];
int iData[1];
numData = 1
iPtr=>iData;
/* get the input data - tag? E? eyp? */ err = OPS_GetIntInput(numData, iPtr)
numData = 2
int numData = 1;
OPS_GetIntInput(&numData, iData); dPtr=>dData;
err = OPS_GetDoubleInput(numData, dPtr)
numData = 2;
OPS_GetDoubleInput(&numData, dData);
c Allocate the element state
matObj%tag = idata(1)
/* Allocate the element state */
matObj%nparam = 2
thisObj->tag = iData[0];
matObj%nstate = 2
thisObj->nParam = 2; /* E, eyp */
What Follows are the Steps
required to Build
[Link] on a
Windows Machine with
Visual Studio 2010 Installed
We Will Build the
[Link]
Source code and example are in
/DEVELOPER/material/cpp
Create Project
1)File>New>Project
4.SelectWin32Project
3.SelectWin32
[Link]
[Link]
[Link]
SelectApplicaAonSeGngs
[Link]
[Link]
[Link]
Add Files To Project
1.
2.
3.
4.
5.
RightClickonSourceFiles
SelectAddExisAng
NavigatetoDEVELOPER/material/cppdirectory
[Link]
SelectAdd
1)SelectBuild>SoluAon
1.
2.
3.
RightClickonElasAcPPcppProject
SelectProperAes
SelectC/C++
IT FAILS!
[Link]
[Link]
1.
2.
SelectC/C++
SelectGeneral
1)SelectBuild>SoluAon
IT FAILS!
1.
2.
3.
4.
5.
RightClickonSourceFiles
SelectAddExisAng
NavigatetoDEVELOPER/coredirectory
[Link]
SelectAdd
[Link]
Copy [Link] from
location into current directory
Now run Example
What Follows are the Steps
required to Build
[Link] on a Linux
Machine with gcc installed
NOTE: We Will use
NEEShub
for this demonstration
([Link]
3 Simple Steps!
1. Download code
svn co svn://[Link]/usr/local/svn/OpenSees/trunk/OpenSees/Developer Developer
2. cd DEVELOPER/material/cpp
3. type make
if you type ls you should see the .so and you
Can test it using >OpenSees [Link]
NOTE: mac users must have xcode installed and
must open DEVELOPER/[Link] and switch
comments on lines 1 and 2 before step 3.
1: Download Source Code
2/3: cd to Directory & Type make
Add the current directory
to your LD_LIBRARY_PATH
Run It