| /trunk/DEVELOPER/coreList |
|---|
| New file |
| 0,0 → 1,149 |
| SRC/analysis/analysis/Analysis.cpp |
| SRC/analysis/analysis/Analysis.h |
| SRC/analysis/model/AnalysisModel.cpp |
| SRC/analysis/model/AnalysisModel.h |
| SRC/tagged/storage/ArrayOfTaggedObjects.cpp |
| SRC/tagged/storage/ArrayOfTaggedObjects.h |
| SRC/tagged/storage/ArrayOfTaggedObjectsIter.cpp |
| SRC/tagged/storage/ArrayOfTaggedObjectsIter.h |
| SRC/handler/BinaryFileStream.cpp |
| SRC/handler/BinaryFileStream.h |
| SRC/actor/channel/Channel.cpp |
| SRC/actor/channel/Channel.h |
| SRC/renderer/ColorMap.h |
| SRC/analysis/handler/ConstraintHandler.cpp |
| SRC/analysis/handler/ConstraintHandler.h |
| SRC/analysis/dof_grp/DOF_Group.cpp |
| SRC/analysis/dof_grp/DOF_Group.h |
| SRC/analysis/model/DOF_GrpIter.cpp |
| SRC/analysis/model/DOF_GrpIter.h |
| SRC/handler/DataFileStream.cpp |
| SRC/handler/DataFileStream.h |
| SRC/domain/domain/Domain.cpp |
| SRC/domain/domain/Domain.h |
| SRC/domain/component/DomainComponent.cpp |
| SRC/domain/component/DomainComponent.h |
| SRC/analysis/analysis/DomainDecompositionAnalysis.h |
| SRC/element/Element.cpp |
| SRC/element/Element.h |
| SRC/domain/domain/ElementIter.h |
| SRC/recorder/response/ElementResponse.cpp |
| SRC/recorder/response/ElementResponse.h |
| SRC/element/ElementalLoad.cpp |
| SRC/element/ElementalLoad.h |
| SRC/domain/load/ElementalLoadIter.cpp |
| SRC/domain/load/ElementalLoadIter.h |
| SRC/analysis/algorithm/equiSolnAlgo/EquiSolnAlgo.h |
| SRC/actor/objectBroker/FEM_ObjectBroker.cpp |
| SRC/actor/objectBroker/FEM_ObjectBroker.h |
| SRC/database/FE_Datastore.cpp |
| SRC/database/FE_Datastore.h |
| SRC/analysis/model/FE_EleIter.cpp |
| SRC/analysis/model/FE_EleIter.h |
| SRC/analysis/fe_ele/FE_Element.cpp |
| SRC/analysis/fe_ele/FE_Element.h |
| SRC//material/section/fiber/Fiber.h |
| SRC/recorder/response/FiberResponse.cpp |
| SRC/recorder/response/FiberResponse.h |
| SRC/G3Globals.h |
| SRC/graph/graph/Graph.cpp |
| SRC/graph/graph/Graph.h |
| SRC/matrix/ID.cpp |
| SRC/matrix/ID.h |
| SRC/analysis/integrator/IncrementalIntegrator.cpp |
| SRC/analysis/integrator/IncrementalIntegrator.h |
| SRC/element/Information.cpp |
| SRC/element/Information.h |
| SRC/analysis/integrator/Integrator.cpp |
| SRC/analysis/integrator/Integrator.h |
| SRC/system_of_eqn/linearSOE/LinearSOE.cpp |
| SRC/system_of_eqn/linearSOE/LinearSOE.h |
| SRC/system_of_eqn/linearSOE/LinearSOESolver.cpp |
| SRC/system_of_eqn/linearSOE/LinearSOESolver.h |
| SRC/domain/load/Load.cpp |
| SRC/domain/load/Load.h |
| SRC/domain/pattern/LoadPattern.h |
| SRC/domain/pattern/LoadPatternIter.cpp |
| SRC/domain/pattern/LoadPatternIter.h |
| SRC/domain/constraints/MP_Constraint.h |
| SRC/domain/constraints/MP_Constraint.cpp |
| SRC/domain/domain/MP_ConstraintIter.h |
| SRC/tagged/storage/MapOfTaggedObjects.cpp |
| SRC/tagged/storage/MapOfTaggedObjects.h |
| SRC/tagged/storage/MapOfTaggedObjectsIter.cpp |
| SRC/tagged/storage/MapOfTaggedObjectsIter.h |
| SRC/material/Material.cpp |
| SRC/material/Material.h |
| SRC/recorder/response/MaterialResponse.cpp |
| SRC/recorder/response/MaterialResponse.h |
| SRC/matrix/Matrix.cpp |
| SRC/matrix/Matrix.h |
| SRC/domain/region/MeshRegion.h |
| SRC/actor/message/Message.cpp |
| SRC/actor/message/Message.h |
| SRC/actor/actor/MovableObject.cpp |
| SRC/actor/actor/MovableObject.h |
| SRC/domain/node/NodalLoad.cpp |
| SRC/domain/node/NodalLoad.h |
| SRC/domain/load/NodalLoadIter.cpp |
| SRC/domain/load/NodalLoadIter.h |
| SRC/domain/node/Node.cpp |
| SRC/domain/node/Node.h |
| SRC/domain/domain/NodeIter.h |
| SRC/OPS_Globals.h |
| SRC/handler/OPS_Stream.cpp |
| SRC/handler/OPS_Stream.h |
| SRC/actor/objectBroker/ObjectBroker.cpp |
| SRC/actor/objectBroker/ObjectBroker.h |
| SRC/domain/component/Parameter.cpp |
| SRC/domain/component/Parameter.h |
| SRC/domain/domain/ParameterIter.h |
| SRC/recorder/Recorder.cpp |
| SRC/recorder/Recorder.h |
| SRC/renderer/Renderer.cpp |
| SRC/renderer/Renderer.h |
| SRC/recorder/response/Response.cpp |
| SRC/recorder/response/Response.h |
| SRC/domain/constraints/SP_Constraint.h |
| SRC/domain/constraints/SP_Constraint.cpp |
| SRC/domain/domain/SP_ConstraintIter.h |
| SRC/domain/domain/single/SingleDomAllSP_Iter.cpp |
| SRC/domain/domain//single/SingleDomAllSP_Iter.h |
| SRC/domain/domain/single/SingleDomEleIter.cpp |
| SRC/domain/domain/single/SingleDomEleIter.h |
| SRC/domain/domain/single/SingleDomMP_Iter.cpp |
| SRC/domain/domain/single/SingleDomMP_Iter.h |
| SRC/domain/domain/single/SingleDomNodIter.cpp |
| SRC/domain/domain/single/SingleDomNodIter.h |
| SRC/domain/domain/single/SingleDomParamIter.cpp |
| SRC/domain/domain/single/SingleDomParamIter.h |
| SRC/domain/domain/single/SingleDomSP_Iter.cpp |
| SRC/domain/domain/single/SingleDomSP_Iter.h |
| SRC/analysis/algorithm/SolutionAlgorithm.cpp |
| SRC/analysis/algorithm/SolutionAlgorithm.h |
| SRC/handler/StandardStream.cpp |
| SRC/handler/StandardStream.h |
| SRC/analysis/integrator/StaticIntegrator.cpp |
| SRC/analysis/integrator/StaticIntegrator.h |
| SRC/domain/subdomain/Subdomain.cpp |
| SRC/domain/subdomain/Subdomain.h |
| SRC/domain/subdomain/SubdomainNodIter.cpp |
| SRC/domain/subdomain/SubdomainNodIter.h |
| SRC/tagged/TaggedObject.cpp |
| SRC/tagged/TaggedObject.h |
| SRC/tagged/storage/TaggedObjectIter.h |
| SRC/tagged/storage/TaggedObjectStorage.h |
| SRC/analysis/integrator/TransientIntegrator.cpp |
| SRC/analysis/integrator/TransientIntegrator.h |
| SRC/material/uniaxial/UniaxialMaterial.cpp |
| SRC/material/uniaxial/UniaxialMaterial.h |
| SRC/matrix/Vector.cpp |
| SRC/matrix/Vector.h |
| SRC/graph/graph/Vertex.cpp |
| SRC/graph/graph/Vertex.h |
| SRC/graph/graph/VertexIter.cpp |
| SRC/graph/graph/VertexIter.h |
| SRC/classTags.h |
| SRC/api/elementAPI.h |
| SRC/coordTransformation/CrdTransf.h |
| PACKAGES/win32/win32Functions.cpp |
| /trunk/DEVELOPER/system/SkylineSPD.cpp |
|---|
| New file |
| 0,0 → 1,622 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision$ |
| // $Date$ |
| // $Source$ |
| // Written: fmk |
| // Revision: A |
| // |
| // Description: This file contains the implementation for ProfileSPDLinSOE |
| #include "SkylineSPD.h" |
| #include <Matrix.h> |
| #include <Graph.h> |
| #include <Vertex.h> |
| #include <VertexIter.h> |
| #include <math.h> |
| #include <Channel.h> |
| #include <FEM_ObjectBroker.h> |
| #include <elementAPI.h> |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export extern "C" _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| // init function called once on first loading of package |
| OPS_Export void |
| localInit() |
| { |
| OPS_Error("SkylineSPD - Written by fmk UC Berkeley Copyright 2009\n", 1); |
| } |
| OPS_Export void * |
| OPS_SkylineSPD() |
| { |
| // |
| // parse args |
| // |
| // none to parse for this soe |
| // |
| // create the SOE |
| // |
| LinearSOE *theSOE = new SkylineSPD(); |
| // return it |
| return theSOE; |
| } |
| // |
| SkylineSPD::SkylineSPD() |
| :LinearSOE(0), |
| size(0), profileSize(0), A(0), B(0), X(0), vectX(0), vectB(0), |
| iDiagLoc(0), Asize(0), Bsize(0), isAfactored(false), isAcondensed(false), |
| numInt(0) , minDiagTol(1e-18), RowTop(0), topRowPtr(0), invD(0) |
| { |
| } |
| SkylineSPD::~SkylineSPD() |
| { |
| if (A != 0) delete [] A; |
| if (B != 0) delete [] B; |
| if (X != 0) delete [] X; |
| if (iDiagLoc != 0) delete [] iDiagLoc; |
| if (vectX != 0) delete vectX; |
| if (vectB != 0) delete vectB; |
| if (RowTop != 0) delete [] RowTop; |
| if (topRowPtr != 0) free((void*)topRowPtr); |
| if (invD != 0) delete [] invD; |
| } |
| int |
| SkylineSPD::getNumEqn(void) const |
| { |
| return size; |
| } |
| int |
| SkylineSPD::setSize(Graph &theGraph) |
| { |
| int oldSize = size; |
| int result = 0; |
| size = theGraph.getNumVertex(); |
| // check we have enough space in iDiagLoc and iLastCol |
| // if not delete old and create new |
| if (size > Bsize) { |
| if (iDiagLoc != 0) delete [] iDiagLoc; |
| if (RowTop != 0) delete [] RowTop; |
| if (topRowPtr != 0) delete [] topRowPtr; |
| if (invD != 0) delete [] invD; |
| // if (topRowPtr != 0) free((void *)topRowPtr); |
| iDiagLoc = new int[size]; |
| RowTop = new int[size]; |
| invD = new double[size]; |
| // topRowPtr = new double *[size] ; |
| topRowPtr = (double **)malloc(size *sizeof(double *)); |
| if (iDiagLoc == 0 || RowTop == 0 || topRowPtr == 0 || invD == 0) { |
| opserr << "WARNING SkylineSPD::setSize() : "; |
| opserr << " - ran out of memory for iDiagLoc\n"; |
| size = 0; Asize = 0; |
| result = -1; |
| } |
| } |
| // zero out iDiagLoc |
| for (int i=0; i<size; i++) { |
| iDiagLoc[i] = 0; |
| } |
| // now we go through the vertices to find the height of each col and |
| // width of each row from the connectivity information. |
| Vertex *vertexPtr; |
| VertexIter &theVertices = theGraph.getVertices(); |
| while ((vertexPtr = theVertices()) != 0) { |
| int vertexNum = vertexPtr->getTag(); |
| const ID &theAdjacency = vertexPtr->getAdjacency(); |
| int iiDiagLoc = iDiagLoc[vertexNum]; |
| int *iiDiagLocPtr = &(iDiagLoc[vertexNum]); |
| for (int i=0; i<theAdjacency.Size(); i++) { |
| int otherNum = theAdjacency(i); |
| int diff = vertexNum-otherNum; |
| if (diff > 0) { |
| if (iiDiagLoc < diff) { |
| iiDiagLoc = diff; |
| *iiDiagLocPtr = diff; |
| } |
| } |
| } |
| } |
| // now go through iDiagLoc, adding 1 for the diagonal element |
| // and then adding previous entry to give current location. |
| if (iDiagLoc != 0) |
| iDiagLoc[0] = 1; // NOTE FORTRAN ARRAY LOCATION |
| for (int j=1; j<size; j++) |
| iDiagLoc[j] = iDiagLoc[j] + 1 + iDiagLoc[j-1]; |
| if (iDiagLoc != 0) |
| profileSize = iDiagLoc[size-1]; |
| // check if we need more space to hold A |
| // if so then go get it |
| if (profileSize > Asize) { |
| // delete old space |
| if (A != 0) |
| delete [] A; |
| // get new space |
| A = new double[profileSize]; |
| if (A == 0) { |
| opserr << "SkylineSPD::SkylineSPD :"; |
| opserr << " ran out of memory for A (size,Profile) ("; |
| opserr << size <<", " << profileSize << ") \n"; |
| size = 0; Asize = 0; profileSize = 0; |
| result = -1; |
| } |
| else |
| Asize = profileSize; |
| } |
| // zero the matrix |
| for (int k=0; k<profileSize; k++) |
| A[k] = 0; |
| isAfactored = false; |
| isAcondensed = false; |
| if (size > Bsize) { // we have to get another space for A |
| // delete the old |
| if (B != 0) delete [] B; |
| if (X != 0) delete [] X; |
| // create the new |
| B = new double[size]; |
| X = new double[size]; |
| if (B == 0 || X == 0 ) { |
| opserr << "SkylineSPD::SkylineSPD :"; |
| opserr << " ran out of memory for vectors (size) ("; |
| opserr << size << ") \n"; |
| size = 0; Bsize = 0; |
| result = -1; |
| } |
| } |
| // zero the vectors |
| for (int l=0; l<size; l++) { |
| B[l] = 0; |
| X[l] = 0; |
| } |
| if (size != oldSize) { |
| if (vectX != 0) |
| delete vectX; |
| if (vectB != 0) |
| delete vectB; |
| vectX = new Vector(X,size); |
| vectB = new Vector(B,size); |
| if (size > Bsize) |
| Bsize = size; |
| } |
| RowTop[0] = 0; |
| topRowPtr[0] = A; |
| for (int j=1; j<size; j++) { |
| int icolsz = iDiagLoc[j] - iDiagLoc[j-1]; |
| RowTop[j] = j - icolsz + 1; |
| topRowPtr[j] = &A[iDiagLoc[j-1]]; // FORTRAN array indexing in iDiagLoc |
| } |
| return result; |
| } |
| int |
| SkylineSPD::addA(const Matrix &m, const ID &id, double fact) |
| { |
| // check for a quick return |
| if (fact == 0.0) return 0; |
| // check that m and id are of similar size |
| int idSize = id.Size(); |
| if (idSize != m.noRows() && idSize != m.noCols()) { |
| opserr << "SkylineSPD::addA() - Matrix and ID not of similar sizes\n"; |
| return -1; |
| } |
| if (fact == 1.0) { // do not need to multiply |
| for (int i=0; i<idSize; i++) { |
| int col = id(i); |
| if (col < size && col >= 0) { |
| double *coliiPtr = &A[iDiagLoc[col] -1]; // -1 as fortran indexing |
| int minColRow; |
| if (col == 0) |
| minColRow = 0; |
| else |
| minColRow = col - (iDiagLoc[col] - iDiagLoc[col-1]) +1; |
| for (int j=0; j<idSize; j++) { |
| int row = id(j); |
| if (row <size && row >= 0 && |
| row <= col && row >= minColRow) { |
| // we only add upper and inside profile |
| double *APtr = coliiPtr + (row-col); |
| *APtr += m(j,i); |
| } |
| } // for j |
| } |
| } // for i |
| } else { |
| for (int i=0; i<idSize; i++) { |
| int col = id(i); |
| if (col < size && col >= 0) { |
| double *coliiPtr = &A[iDiagLoc[col] -1]; // -1 as fortran indexing |
| int minColRow; |
| if (col == 0) |
| minColRow = 0; |
| else |
| minColRow = col - (iDiagLoc[col] - iDiagLoc[col-1]) +1; |
| for (int j=0; j<idSize; j++) { |
| int row = id(j); |
| if (row <size && row >= 0 && |
| row <= col && row >= minColRow) { |
| // we only add upper and inside profile |
| double *APtr = coliiPtr + (row-col); |
| *APtr += m(j,i) * fact; |
| } |
| } // for j |
| } |
| } // for i |
| } |
| return 0; |
| } |
| int |
| SkylineSPD::addB(const Vector &v, const ID &id, double fact) |
| { |
| // check for a quick return |
| if (fact == 0.0) return 0; |
| // check that m and id are of similar size |
| int idSize = id.Size(); |
| if (idSize != v.Size() ) { |
| opserr << "SkylineSPD::addB() -"; |
| opserr << " Vector and ID not of similar sizes\n"; |
| return -1; |
| } |
| if (fact == 1.0) { // do not need to multiply if fact == 1.0 |
| for (int i=0; i<id.Size(); i++) { |
| int pos = id(i); |
| if (pos <size && pos >= 0) |
| B[pos] += v(i); |
| } |
| } else if (fact == -1.0) { // do not need to multiply if fact == -1.0 |
| for (int i=0; i<id.Size(); i++) { |
| int pos = id(i); |
| if (pos <size && pos >= 0) |
| B[pos] -= v(i); |
| } |
| } else { |
| for (int i=0; i<id.Size(); i++) { |
| int pos = id(i); |
| if (pos <size && pos >= 0) |
| B[pos] += v(i) * fact; |
| } |
| } |
| return 0; |
| } |
| int |
| SkylineSPD::setB(const Vector &v, double fact) |
| { |
| // check for a quick return |
| if (fact == 0.0) return 0; |
| if (v.Size() != size) { |
| opserr << "WARNING BandGenLinSOE::setB() -"; |
| opserr << " incomptable sizes " << size << " and " << v.Size() << endln; |
| return -1; |
| } |
| if (fact == 1.0) { // do not need to multiply if fact == 1.0 |
| for (int i=0; i<size; i++) { |
| B[i] = v(i); |
| } |
| } else if (fact == -1.0) { |
| for (int i=0; i<size; i++) { |
| B[i] = -v(i); |
| } |
| } else { |
| for (int i=0; i<size; i++) { |
| B[i] = v(i) * fact; |
| } |
| } |
| return 0; |
| } |
| void |
| SkylineSPD::zeroA(void) |
| { |
| double *Aptr = A; |
| for (int i=0; i<Asize; i++) |
| *Aptr++ = 0; |
| isAfactored = false; |
| } |
| void |
| SkylineSPD::zeroB(void) |
| { |
| double *Bptr = B; |
| for (int i=0; i<size; i++) |
| *Bptr++ = 0; |
| } |
| void |
| SkylineSPD::setX(int loc, double value) |
| { |
| if (loc < size && loc >=0) |
| X[loc] = value; |
| } |
| void |
| SkylineSPD::setX(const Vector &x) |
| { |
| if (x.Size() == size && vectX != 0) |
| *vectX = x; |
| } |
| const Vector & |
| SkylineSPD::getX(void) |
| { |
| if (vectX == 0) { |
| opserr << "FATAL SkylineSPD::getX - vectX == 0"; |
| exit(-1); |
| } |
| return *vectX; |
| } |
| const Vector & |
| SkylineSPD::getB(void) |
| { |
| if (vectB == 0) { |
| opserr << "FATAL SkylineSPD::getB - vectB == 0"; |
| exit(-1); |
| } |
| return *vectB; |
| } |
| double |
| SkylineSPD::normRHS(void) |
| { |
| double norm =0.0; |
| for (int i=0; i<size; i++) { |
| double Yi = B[i]; |
| norm += Yi*Yi; |
| } |
| return sqrt(norm); |
| } |
| int |
| SkylineSPD::sendSelf(int cTag, Channel &theChannel) |
| { |
| return 0; |
| } |
| int |
| SkylineSPD::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker) |
| { |
| return 0; |
| } |
| int |
| SkylineSPD::solve(void) |
| { |
| if (size == 0) |
| return 0; |
| // copy B into X |
| for (int ii=0; ii<size; ii++) |
| X[ii] = B[ii]; |
| if (isAfactored == false) { |
| // FACTOR & SOLVE |
| double *ajiPtr, *akjPtr, *akiPtr, *bjPtr; |
| // if the matrix has not been factored already factor it into U^t D U |
| // storing D^-1 in invD as we go |
| double a00 = A[0]; |
| if (a00 <= 0.0) { |
| opserr << "ProfileSPDLinDirectSolver::solve() - "; |
| opserr << " aii < 0 (i, aii): (0,0)\n"; |
| return(-2); |
| } |
| invD[0] = 1.0/A[0]; |
| // for every col across |
| for (int i=1; i<size; i++) { |
| int rowitop = RowTop[i]; |
| ajiPtr = topRowPtr[i]; |
| for (int j=rowitop; j<i; j++) { |
| double tmp = *ajiPtr; |
| int rowjtop = RowTop[j]; |
| if (rowitop > rowjtop) { |
| akjPtr = topRowPtr[j] + (rowitop-rowjtop); |
| akiPtr = topRowPtr[i]; |
| for (int k=rowitop; k<j; k++) |
| tmp -= *akjPtr++ * *akiPtr++ ; |
| *ajiPtr++ = tmp; |
| } |
| else { |
| akjPtr = topRowPtr[j]; |
| akiPtr = topRowPtr[i] + (rowjtop-rowitop); |
| for (int k=rowjtop; k<j; k++) |
| tmp -= *akjPtr++ * *akiPtr++ ; |
| *ajiPtr++ = tmp; |
| } |
| } |
| /* now form i'th col of [U] and determine [dii] */ |
| double aii = A[iDiagLoc[i] -1]; // FORTRAN ARRAY INDEXING |
| ajiPtr = topRowPtr[i]; |
| double *bjPtr = &X[rowitop]; |
| double tmp = 0; |
| for (int jj=rowitop; jj<i; jj++) { |
| double aji = *ajiPtr; |
| double lij = aji * invD[jj]; |
| tmp -= lij * *bjPtr++; |
| *ajiPtr++ = lij; |
| aii = aii - lij*aji; |
| } |
| // check that the diag > the tolerance specified |
| if (aii == 0.0) { |
| opserr << "ProfileSPDLinDirectSolver::solve() - "; |
| opserr << " aii < 0 (i, aii): (" << i << ", " << aii << ")\n"; |
| return(-2); |
| } |
| if (fabs(aii) <= minDiagTol) { |
| opserr << "ProfileSPDLinDirectSolver::solve() - "; |
| opserr << " aii < minDiagTol (i, aii): (" << i; |
| opserr << ", " << aii << ")\n"; |
| return(-2); |
| } |
| invD[i] = 1.0/aii; |
| X[i] += tmp; |
| } |
| isAfactored = true; |
| numInt = 0; |
| // divide by diag term |
| bjPtr = X; |
| double *aiiPtr = invD; |
| for (int j=0; j<size; j++) |
| *bjPtr++ = *aiiPtr++ * X[j]; |
| // now do the back substitution storing result in X |
| for (int k=(size-1); k>0; k--) { |
| int rowktop = RowTop[k]; |
| double bk = X[k]; |
| double *ajiPtr = topRowPtr[k]; |
| for (int j=rowktop; j<k; j++) |
| X[j] -= *ajiPtr++ * bk; |
| } |
| } |
| else { |
| // JUST DO SOLVE |
| // do forward substitution |
| for (int i=1; i<size; i++) { |
| int rowitop = RowTop[i]; |
| double *ajiPtr = topRowPtr[i]; |
| double *bjPtr = &X[rowitop]; |
| double tmp = 0; |
| for (int j=rowitop; j<i; j++) |
| tmp -= *ajiPtr++ * *bjPtr++; |
| X[i] += tmp; |
| } |
| // divide by diag term |
| double *bjPtr = X; |
| double *aiiPtr = invD; |
| for (int j=0; j<size; j++) |
| *bjPtr++ = *aiiPtr++ * X[j]; |
| // now do the back substitution storing result in X |
| for (int k=(size-1); k>0; k--) { |
| int rowktop = RowTop[k]; |
| double bk = X[k]; |
| double *ajiPtr = topRowPtr[k]; |
| for (int j=rowktop; j<k; j++) |
| X[j] -= *ajiPtr++ * bk; |
| } |
| } |
| return 0; |
| } |
| /trunk/DEVELOPER/system/SkylineSPD.h |
|---|
| New file |
| 0,0 → 1,90 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision$ |
| // $Date$ |
| // $Source$ |
| // Written: fmk |
| // Revision: A |
| // |
| // Description: This file contains the class definition for MyProfileSPD_SOE |
| // SkylineSPD is a subclass of LinearSOE. It stores symmetric positive |
| // definite matrix in skyline and solves directly using a cholesky LL^t |
| // decomposition. |
| // What: "@(#) SkylineSPD.h, revA" |
| #ifndef SkylineSPD_h |
| #define SkylineSPD_h |
| #include <LinearSOE.h> |
| #include <Vector.h> |
| class SkylineSPD : public LinearSOE |
| { |
| public: |
| SkylineSPD(); |
| ~SkylineSPD(); |
| int getNumEqn(void) const; |
| int setSize(Graph &theGraph); |
| int addA(const Matrix &, const ID &, double fact = 1.0); |
| int addB(const Vector &, const ID &, double fact = 1.0); |
| int setB(const Vector &, double fact = 1.0); |
| void zeroA(void); |
| void zeroB(void); |
| void setX(int loc, double value); |
| void setX(const Vector &x); |
| const Vector &getX(void); |
| const Vector &getB(void); |
| double normRHS(void); |
| int solve(void); |
| int sendSelf(int commitTag, Channel &theChannel); |
| int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker); |
| protected: |
| private: |
| int size, profileSize; |
| double *A, *B, *X; |
| Vector *vectX; |
| Vector *vectB; |
| int *iDiagLoc; |
| int Asize, Bsize; |
| bool isAfactored, isAcondensed; |
| int numInt; |
| double minDiagTol; |
| int *RowTop; |
| double **topRowPtr, *invD; |
| }; |
| #endif |
| /trunk/DEVELOPER/system/Makefile |
|---|
| New file |
| 0,0 → 1,42 |
| include ../Makefile.def |
| OS_FLAG ?= -D_LINUX |
| ifeq ($(OS_FLAG), -D_LINUX) |
| LIBRARY = SkylineSPD.so |
| all: $(LIBRARY) |
| OBJS = SkylineSPD.o |
| $(LIBRARY): $(OBJS) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY) $(OBJS) -lc |
| SkylineSPD.o: SkylineSPD.cpp |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall SkylineSPD.cpp |
| else |
| LIBRARY = SkylineSPD.dylib |
| all: |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace SkylineSPD.cpp $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY) |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/system/example1.tcl |
|---|
| New file |
| 0,0 → 1,65 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.2 $ |
| # $Date: 2008/12/09 22:43:40 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewElement/cpp/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial Elastic 1 3000 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element truss 1 1 4 10.0 1 |
| element truss 2 2 4 5.0 1 |
| element truss 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| #system ProfileSPD |
| system SkylineSPD |
| constraints Plain |
| integrator LoadControl 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| system ProfileSPD |
| analyze 1 |
| system SkylineSPD |
| analyze 1 |
| exit |
| /trunk/DEVELOPER/recorder/SumElementForcesRecorder.cpp |
|---|
| New file |
| 0,0 → 1,427 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision$ |
| // $Date$ |
| // $Source$ |
| // Written: fmk |
| // Created: 05/10 |
| // |
| // Description: This file contains the class implementatation of SumElementForcesRecorder. |
| // |
| // What: "@(#) SumElementForcesRecorder.C, revA" |
| #include "SumElementForcesRecorder.h" |
| #include <elementAPI.h> |
| #include <Domain.h> |
| #include <Element.h> |
| #include <ElementIter.h> |
| #include <Matrix.h> |
| #include <Vector.h> |
| #include <ID.h> |
| #include <string.h> |
| #include <Response.h> |
| #include <Message.h> |
| #include <Channel.h> |
| #include <FEM_ObjectBroker.h> |
| #include <StandardStream.h> |
| #include <BinaryFileStream.h> |
| #include <DataFileStream.h> |
| #include <elementAPI.h> |
| SumElementForcesRecorder::SumElementForcesRecorder() |
| :Recorder(-1), |
| numEle(0), theElements(0), eleID(0), |
| theDomain(0), theOutput(0), |
| echoTimeFlag(true), data(0) |
| { |
| } |
| SumElementForcesRecorder::SumElementForcesRecorder(const ID ele, |
| bool echoTime, |
| OPS_Stream *theoutput) |
| :Recorder(-1), |
| numEle(0), theElements(0), eleID(ele), |
| theDomain(0), theOutput(theoutput), |
| echoTimeFlag(echoTime), data(0) |
| { |
| } |
| SumElementForcesRecorder::~SumElementForcesRecorder() |
| { |
| if (theElements != 0) |
| delete [] theElements; |
| if (data != 0) |
| delete data; |
| if (theOutput != 0) |
| delete theOutput; |
| } |
| int |
| SumElementForcesRecorder::record(int commitTag, double timeStamp) |
| { |
| // check for initialization |
| if (data == 0) { |
| opserr << "SumElementForcesRecorder::record() - setDomain() has not been called\n"; |
| return -1; |
| } |
| // zero the data vector |
| data->Zero(); |
| int forceSize = data->Size(); |
| int startLoc = 0; |
| // write the time if echTimeFlag set |
| if (echoTimeFlag == true) { |
| (*data)(0) = timeStamp; |
| forceSize -= 1; |
| startLoc = 1; |
| } |
| // |
| // for each element that has been added to theElements add force contribution |
| // |
| for (int i=0; i< numEle; i++) { |
| if (theElements[i] != 0) { |
| int loc = startLoc; |
| const Vector &force = theElements[i]->getResistingForce(); |
| int forceSize = force.Size(); |
| for (int j=0; j<forceSize; j++, loc++) |
| (*data)(loc) += force(j); |
| } |
| } |
| // |
| // send the response vector to the output handler for o/p |
| // |
| theOutput->write(*data); |
| // succesfull completion - return 0 |
| return 0; |
| } |
| int |
| SumElementForcesRecorder::restart(void) |
| { |
| return 0; |
| } |
| int |
| SumElementForcesRecorder::domainChanged(void) |
| { |
| if (theDomain != 0) |
| this->setDomain(*theDomain); |
| } |
| int |
| SumElementForcesRecorder::setDomain(Domain &theDom) |
| { |
| theDomain = &theDom; |
| // delete old data |
| if (theElements != 0) |
| delete [] theElements; |
| if (data != 0) |
| delete data; |
| // set numEle |
| numEle = eleID.Size(); |
| if (numEle == 0) { |
| opserr << "WARNING SumElementForcesRecorder::initialize() - no elements tags passed in input!\n"; |
| return 0; |
| } |
| // create theElements, an array of pointers to elements |
| theElements = new Element *[numEle]; |
| if (theElements == 0) { |
| opserr << "WARNING SumElementForcesRecorder::initialize() - out of memory\n"; |
| numEle = 0; // set numEle = 0, in case record() still called |
| return -1; |
| } |
| // zero theElements array |
| for (int k=0; k<numEle; k++) |
| theElements[k] = 0; |
| // |
| // loop over the list of elements, |
| // if element exists add it's pointer o the array |
| // get its resisting force, check size to determine compatable with others |
| // |
| int sizeArray = -1; |
| for (int i=0; i<numEle; i++) { |
| int eleTag = eleID(i); |
| Element *theEle = theDomain->getElement(eleTag); |
| if (theEle != 0) { |
| const Vector &force = theEle->getResistingForce(); |
| int forceSize = force.Size(); |
| if (sizeArray == -1) { |
| sizeArray = forceSize; |
| theElements[i] = theEle; |
| } else if (sizeArray != forceSize) { |
| opserr << "WARNING: forces mismatch - element: " << eleTag << " will not be included\n"; |
| } else { |
| theElements[i] = theEle; |
| } |
| } |
| } |
| // if echTimeFlag is set, add room for the time to be output |
| if (echoTimeFlag == true) |
| sizeArray++; |
| // create the vector to hold the data |
| data = new Vector(sizeArray); |
| if (data == 0 || data->Size() != sizeArray) { |
| opserr << "SumElementForcesRecorder::initialize() - out of memory\n"; |
| delete [] theElements; |
| theElements = 0; |
| numEle = 0; |
| } |
| return 0; |
| } |
| static char myClassType[] = {"SumElementForcesRecorder"}; |
| const char * |
| SumElementForcesRecorder::getClassType(void) const |
| { |
| return myClassType; |
| } |
| int |
| SumElementForcesRecorder::sendSelf(int commitTag, Channel &theChannel) |
| { |
| // send in an ID (integar array) to the receiving object the following: |
| // recorder tag |
| // size of eleID |
| // class tag of handler |
| // echoTimeFlag |
| static ID idData(5); |
| idData(0) = this->getTag();; |
| idData(1) = eleID.Size(); |
| idData(2) = theOutput->getClassTag(); |
| if (echoTimeFlag == true) |
| idData(3) = 1; |
| else |
| idData(3) = 0; |
| if (theChannel.sendID(0, commitTag, idData) < 0) { |
| opserr << "SumElementForcesRecorder::recvSelf() - failed to recv idData\n"; |
| return -1; |
| } |
| // send eleID to receiving object |
| if (theChannel.sendID(0, commitTag, eleID) < 0) { |
| opserr << "SumElementForcesRecorder::sendSelf() - failed to send idData\n"; |
| return -1; |
| } |
| // send theOutput to receiving object |
| if (theOutput->sendSelf(commitTag, theChannel) < 0) { |
| opserr << "SumElementForcesRecorder::sendSelf() - failed to send theOutput\n"; |
| return -1; |
| } |
| return 0; |
| } |
| int |
| SumElementForcesRecorder::recvSelf(int commitTag, Channel &theChannel, |
| FEM_ObjectBroker &theBroker) |
| { |
| // receive from the sending object the ID |
| static ID idData(5); |
| if (theChannel.recvID(0, commitTag, idData) < 0) { |
| opserr << "SumElementForcesRecorder::recvSelf() - failed to recv idData\n"; |
| return -1; |
| } |
| // with the data received |
| // setTag |
| // resize the eleID array |
| // set echoTimeFlag |
| // get an outputHandler |
| this->setTag(idData(0)); |
| eleID.resize(idData(1)); |
| idData(2) = theOutput->getClassTag(); |
| if (idData(3) == 0) |
| echoTimeFlag = true; |
| else |
| echoTimeFlag = false; |
| if (theOutput != 0 && theOutput->getClassTag() != idData(4)) |
| delete theOutput; |
| theOutput = theBroker.getPtrNewStream(idData(4)); |
| if (theOutput == 0) { |
| opserr << "SumElementForcesRecorder::recvSelf() - failed to get Output of correct type\n"; |
| return -1; |
| } |
| // receive eleID |
| if (theChannel.recvID(0, commitTag, eleID) < 0) { |
| opserr << "SumElementForcesRecorder::recvSelf() - failed to recv eleID\n"; |
| return -1; |
| } |
| // get theOutput to receive data |
| if (theOutput->recvSelf(commitTag, theChannel, theBroker) < 0) { |
| opserr << "SumElementForcesRecorder::sendSelf() - failed to send theOutput\n"; |
| return -1; |
| } |
| return 0; |
| } |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export extern "C" _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| static int numSumElementForcesREcorder = 0; |
| OPS_Export void * |
| OPS_SumElementForcesRecorder() |
| { |
| Recorder *theRecorder = 0; |
| int numRemainingArgs = OPS_GetNumRemainingInputArgs(); |
| // check for quick return, possibly parallel case |
| if (numRemainingArgs == 0) { |
| Recorder *theRecorder = new SumElementForcesRecorder(); |
| } |
| // |
| // parse args |
| // |
| int numEle = 0, eleTag; |
| ID eleID(0); |
| OPS_Stream *theOutputStream = 0; |
| int outMode = 0; // standard stream |
| bool echoTime = false; |
| bool doneParsingArgs = false; |
| char data[100]; |
| char outputName[200]; |
| char **eleArgs = 0; |
| int numEleArgs = 0; |
| while (numRemainingArgs > 0) { |
| if (OPS_GetString(data,100) < 0) |
| return 0; |
| // output to standard file |
| if (strcmp(data,"-file") == 0) { |
| outMode = 1; |
| if (OPS_GetString(outputName,200) < 0) |
| return 0; |
| numRemainingArgs -= 2; |
| } |
| // output to binary file |
| else if (strcmp(data,"-binary") == 0) { |
| outMode = 2; |
| if (OPS_GetString(outputName,200) < 0) |
| return 0; |
| numRemainingArgs -= 2; |
| } |
| // echo domain time stamp in output |
| else if (strcmp(data,"-time") == 0) { |
| echoTime = true; |
| numRemainingArgs -= 1; |
| } |
| // read the list of elements & place in an ID |
| else if ((strcmp(data,"-ele") == 0) || |
| (strcmp(data,"-eles") == 0) || |
| (strcmp(data,"-element") == 0)) { |
| numRemainingArgs --; |
| int one = 1; |
| while (numRemainingArgs > 0 && OPS_GetIntInput(&one, &eleTag) == 0) { |
| eleID[numEle] = eleTag; |
| numEle++; |
| numRemainingArgs--; |
| } |
| doneParsingArgs = true; |
| } |
| else |
| doneParsingArgs = true; |
| } |
| // |
| // create the output handler |
| // |
| if (outMode == 0) |
| theOutputStream = new StandardStream(); |
| if (outMode == 1) |
| theOutputStream = new DataFileStream(outputName); |
| else if (outMode == 2) |
| theOutputStream = new BinaryFileStream(outputName); |
| // |
| // create the recorder |
| // |
| theRecorder = new SumElementForcesRecorder(eleID, |
| echoTime, |
| theOutputStream); |
| // return it |
| return theRecorder; |
| } |
| /trunk/DEVELOPER/recorder/SumElementForcesRecorder.h |
|---|
| New file |
| 0,0 → 1,81 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision: 1.1 $ |
| // $Date: 2009/03/20 22:44:17 $ |
| // $Source: /usr/local/cvs/OpenSees/PACKAGES/NewRecorder/SumElementForcesRecorder.h,v $ |
| #ifndef SumElementForcesRecorder_h |
| #define SumElementForcesRecorder_h |
| // Written: fmk |
| // Revision: A |
| // |
| // Description: This file contains the class definition for SumElementForcesRecorder. |
| // A SumElementForcesRecorder is used to obtain a response from an element during |
| // the analysis. |
| // |
| // What: "@(#) SumElementForcesRecorder.h, revA" |
| #include <Recorder.h> |
| #include <Information.h> |
| #include <ID.h> |
| class Domain; |
| class Vector; |
| class Matrix; |
| class Element; |
| class Response; |
| class FE_Datastore; |
| class SumElementForcesRecorder: public Recorder |
| { |
| public: |
| SumElementForcesRecorder(); |
| SumElementForcesRecorder(const ID eleID, |
| bool echoTime, |
| OPS_Stream *theOutputHandler); |
| ~SumElementForcesRecorder(); |
| int record(int commitTag, double timeStamp); |
| int restart(void); |
| int domainChanged(void); |
| int setDomain(Domain &theDomain); |
| const char *getClassType(void) const; |
| int sendSelf(int commitTag, Channel &theChannel); |
| int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker); |
| protected: |
| private: |
| int numEle; // the number of elements |
| Element **theElements;// pointer to array of element pointers |
| ID eleID; // ID (integer list) of element tags to record |
| Domain *theDomain; // pointer to domain holding elements |
| OPS_Stream *theOutput;// pointer to output location |
| bool echoTimeFlag; // flag indicating if pseudo time to be printed |
| Vector *data; // Vector (double array) to store sum of element forces |
| }; |
| #endif |
| /trunk/DEVELOPER/recorder/Makefile |
|---|
| New file |
| 0,0 → 1,40 |
| include ../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| LIBRARY = SumElementForcesRecorder.so |
| all: $(LIBRARY) |
| OBJS = SumElementForcesRecorder.o |
| $(LIBRARY): $(OBJS) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY) $(OBJS) -lc |
| SumElementForcesRecorder.o: SumElementForcesRecorder.cpp |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall SumElementForcesRecorder.cpp |
| else |
| LIBRARY = SumElementForcesRecorder.dylib |
| all: |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace SumElementForcesRecorder.cpp $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY) |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/recorder/example1.tcl |
|---|
| New file |
| 0,0 → 1,60 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.1 $ |
| # $Date: 2009/03/20 22:44:17 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewRecorder/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial Elastic 1 3000 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element truss 1 1 4 10.0 1 |
| element truss 2 2 4 5.0 1 |
| element truss 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system ProfileSPD |
| constraints Plain |
| integrator LoadControl 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Element -file a.out -time -ele 1 2 3 forces |
| recorder SumElementForcesRecorder -file b.out -time -ele 1 2 3 |
| #recorder ElementSumRecorder -file b.out -time -ele 1 2 3 forces |
| # perform the analysis |
| analyze 10 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| exit |
| /trunk/DEVELOPER/makeCore.tcl |
|---|
| New file |
| 0,0 → 1,36 |
| set openSeesDir "/Users/fmk/OpenSees" |
| set coreDir "/Users/fmk/OpenSees/DEVELOPER/core" |
| set coreList "/Users/fmk/OpenSees/DEVELOPER/coreList"; |
| # procedure to go through the directories and clean out stuff |
| proc copyIt {openSeesDir fileName coreDir} { |
| set dirFile $openSeesDir |
| set dirFile [file join $dirFile $fileName] |
| file copy -force $dirFile $coreDir |
| } |
| proc makeCore {openSeesDir coreList coreDir} { |
| file mkdir $coreDir |
| set pwd [pwd] |
| set pwd [file join $pwd "core"] |
| set inFileID [open $coreList r] |
| foreach line [split [read $inFileID] \n] { |
| if {[llength $line] == 0} { |
| # Blank line --> do nothing |
| continue |
| } else { |
| copyIt $openSeesDir $line $coreDir |
| } |
| } |
| puts "DONE copying files from $coreList to $coreDir" |
| } |
| makeCore $openSeesDir $coreList $coreDir |
| /trunk/DEVELOPER/element/c/trussC.c |
|---|
| New file |
| 0,0 → 1,239 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| /* |
| ** $Revision: 1.6 $ |
| ** $Date: 2009/01/10 00:33:54 $ |
| ** $Source: /usr/local/cvs/OpenSees/PACKAGES/NewElement/c/trussC.c,v $ |
| ** Written: fmk |
| ** |
| ** Description: This file contains the implementation of truss element. |
| */ |
| #include "elementAPI.h" |
| #include <math.h> |
| #include <stdlib.h> |
| #include <stdio.h> |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| OPS_Export void |
| trussC (eleObj *thisObj, modelState *model, double *tang, double *resid, int *isw, int *errFlag) |
| { |
| double matStrain[1]; |
| double matTang[1]; |
| double matStress[1]; |
| double A, L, K, rho, cs, sn; |
| int i,j; |
| double dx, dy, dLength; |
| double d1[2]; |
| double d2[2]; |
| int nd1, nd2, numDOF, numCrd; |
| int *matData, matTag, matType, n1, n2; |
| double nd1Crd[2]; |
| double nd2Crd[2]; |
| double tran[4]; |
| matObj *theMat; |
| if (*isw == ISW_INIT) { |
| double dData[1]; |
| int iData[4]; |
| /* get the input data - tag? nd1? nd2? A? matTag? */ |
| int numData = 3; |
| OPS_GetIntInput(&numData, iData); |
| numData = 1; |
| OPS_GetDoubleInput(&numData, dData); |
| OPS_GetIntInput(&numData, &iData[3]); |
| /* Allocate the element state */ |
| thisObj->tag = iData[0]; |
| n1 = iData[1]; |
| n2 = iData[2]; |
| matTag = iData[3]; |
| thisObj->nNode = 2; |
| thisObj->nParam = 4; |
| thisObj->nDOF = 4; |
| thisObj->nState = 0; |
| thisObj->nMat = 1; |
| iData[0] = matTag; |
| matData = iData; |
| matType = OPS_UNIAXIAL_MATERIAL_TYPE; |
| OPS_AllocateElement(thisObj, matData, &matType); |
| /* fill in the element state */ |
| thisObj->node[0] = n1; |
| thisObj->node[1] = n2; |
| numCrd = 2; |
| thisObj->param[0] = dData[0]; |
| OPS_GetNodeCrd(&n1, &numCrd, nd1Crd); |
| OPS_GetNodeCrd(&n2, &numCrd, nd2Crd); |
| dx = nd2Crd[0]-nd1Crd[0]; |
| dy = nd2Crd[1]-nd1Crd[1]; |
| L = sqrt(dx*dx + dy*dy); |
| thisObj->param[1] = L; |
| if (L == 0.0) { |
| OPS_Error("Warning - truss element has zero length\n", 1); |
| return; |
| } |
| cs = dx/L; |
| sn = dy/L; |
| thisObj->param[2] = cs; |
| thisObj->param[3] = sn; |
| /* ****************************** |
| placed in AllocateElement |
| matObj *theMat = OPS_GetMaterial(&(iData[3])); |
| if (theMat == 0) { |
| // OPS_Error("Warning - truss element could not find material\n",1); |
| return; |
| } |
| thisObj->mats[0] = theMat; |
| ********************************************** */ |
| *errFlag = 0; |
| } else if (*isw == ISW_COMMIT) { |
| matObj *theMat = thisObj->mats[0]; |
| theMat->matFunctPtr(theMat, model, matStrain, matTang, matStress, isw, errFlag); |
| } else if (*isw == ISW_REVERT) { |
| matObj *theMat = thisObj->mats[0]; |
| theMat->matFunctPtr(theMat, model, matStrain, matTang, matStress, isw, errFlag); |
| } else if (*isw == ISW_REVERT_TO_START) { |
| matObj *theMat = thisObj->mats[0]; |
| theMat->matFunctPtr(theMat, model, matStrain, matTang, matStress, isw, errFlag); |
| } else if (*isw == ISW_FORM_TANG_AND_RESID) { |
| double L = thisObj->param[1]; |
| if (L == 0.0) { |
| // OPS_Error("Warning - truss element has zero length\n", 1); |
| return; |
| } |
| nd1 = thisObj->node[0]; |
| nd2 = thisObj->node[1]; |
| numDOF = 2; |
| OPS_GetNodeDisp(&nd1, &numDOF, d1); |
| OPS_GetNodeDisp(&nd2, &numDOF, d2); |
| A = thisObj->param[0]; |
| cs = thisObj->param[2]; |
| sn = thisObj->param[3]; |
| tran[0] = -cs; |
| tran[1] = -sn; |
| tran[2] = cs; |
| tran[3] = sn; |
| dLength = 0.0; |
| for (i=0; i<2; i++){ |
| dLength -= (d2[i]-d1[i]) * tran[i]; |
| } |
| matStrain[0] = dLength/L; |
| theMat = thisObj->mats[0]; |
| theMat->matFunctPtr(theMat, model, matStrain, matTang, matStress, isw, errFlag); |
| /* ******************* instead of call material funtion |
| int matNum = 0; |
| *errFlag = OPS_InvokeMaterial(thisObj, &matNum, model, matStrain, matStress, matTang, isw); |
| fprintf(stderr,"strain, tang, stress: %e %e %e\n",matStrain[0], matTang[0], matStress[0]); |
| *errFlag = OPS_InvokeMaterialDirectly(theMat, model, matStrain, matStress, matTang, isw); |
| ******************************************************** */ |
| if (*errFlag == 0) { |
| double force = A*matStress[0]; |
| for (i=0; i<4; i++) |
| resid[i] = tran[i]*force; |
| K = A*matTang[0]/L; |
| // tang(j,i) |
| for (i = 0; i<4; i++) |
| for (j=0; j < 4; j++) |
| tang[i+j*4] = K * tran[i]*tran[j]; |
| } |
| } else if (*isw == ISW_FORM_MASS) { |
| L = thisObj->param[1]; |
| if (L == 0.0) { |
| // OPS_Error("Warning - truss element has zero length\n", 1); |
| return; |
| } |
| A = thisObj->param[0]; |
| rho = thisObj->param[4]; |
| for (i=0; i<16; i++) |
| tang[i] = 0.0; |
| if (rho != 0.0) { |
| double massV = rho * A * L/2; |
| tang[0] = massV; |
| tang[1+1*4] = massV; |
| tang[2+2*4] = massV; |
| tang[3+3*4] = massV; |
| } |
| *errFlag = 0; |
| } |
| } |
| OPS_Export void |
| localInit() |
| { |
| OPS_Error("trussC element - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n", 1); |
| } |
| /trunk/DEVELOPER/element/c/Makefile |
|---|
| New file |
| 0,0 → 1,39 |
| include ../../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| all: $(LIBRARY_C) |
| LIBRARY_C = trussC.so |
| OBJS_C = TrussC.o |
| $(LIBRARY_C): $(OBJS_C) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY_CPP) $(OBJS_CPP) -lc |
| trussC.o: trussC.c |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall trussC.c |
| else |
| LIBRARY_C = trussC.dylib |
| all: $(OBJS) |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace trussC.c $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY_C) |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/element/c/example1.tcl |
|---|
| New file |
| 0,0 → 1,62 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.1 $ |
| # $Date: 2008/12/01 23:30:25 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewElement/c/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial Elastic 1 3000 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element trussC 1 1 4 10.0 1 |
| element trussC 2 2 4 5.0 1 |
| element trussC 3 3 4 5.0 1 |
| #element truss 1 1 4 10.0 1 |
| #element truss 2 2 4 5.0 1 |
| #element truss 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system BandSPD |
| constraints Plain |
| integrator LoadControl 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| /trunk/DEVELOPER/element/cpp/Truss2D.cpp |
|---|
| New file |
| 0,0 → 1,553 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // Written: fmk |
| // |
| // What: "@(#) Truss2D.C, revA" |
| // we specify what header files we need |
| #include "Truss2D.h" |
| #include <elementAPI.h> |
| #include <G3Globals.h> |
| #include <Information.h> |
| #include <Domain.h> |
| #include <Node.h> |
| #include <Channel.h> |
| #include <Message.h> |
| #include <FEM_ObjectBroker.h> |
| #include <UniaxialMaterial.h> |
| #include <Renderer.h> |
| #include <ElementResponse.h> |
| #include <math.h> |
| #include <stdlib.h> |
| #include <string.h> |
| // initialise the class wide variables |
| Matrix Truss2D::trussK(4,4); |
| Vector Truss2D::trussR(4); |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export extern "C" _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| static int numMyTruss = 0; |
| OPS_Export void * |
| OPS_Truss2D() |
| { |
| // print out a message about who wrote this element & any copyright info wanted |
| if (numMyTruss == 0) { |
| opserr << "Truss2D element - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n"; |
| numMyTruss++; |
| } |
| Element *theTruss = 0; |
| int numRemainingArgs = OPS_GetNumRemainingInputArgs(); |
| if (numRemainingArgs == 0) { // parallel processing |
| theTruss = new Truss2D(); |
| return theTruss; |
| } |
| if (numRemainingArgs != 5) { |
| opserr << "ERROR - Truss2D not enough args provided, want: element Truss2D tag? iNode? jNode? Area? matTag?\n"; |
| numMyTruss++; |
| } |
| // get the id and end nodes |
| int iData[4]; |
| double dData[1]; |
| int numData; |
| numData = 3; |
| if (OPS_GetIntInput(&numData, iData) != 0) { |
| opserr << "WARNING invalid element data\n"; |
| return 0; |
| } |
| int eleTag = iData[0]; |
| numData = 1; |
| if (OPS_GetDoubleInput(&numData, dData) != 0) { |
| opserr << "WARNING error reading element area for element" << eleTag << endln; |
| return 0; |
| } |
| numData = 1; |
| if (OPS_GetIntInput(&numData, &iData[3]) != 0) { |
| opserr << "WARNING error reading element material tag for element " << eleTag << endln; |
| return 0; |
| } |
| int matID = iData[3]; |
| UniaxialMaterial *theMaterial = OPS_GetUniaxialMaterial(matID); |
| if (theMaterial == 0) { |
| opserr << "WARNING material with tag " << matID << "not found for element " << eleTag << endln; |
| return 0; |
| } |
| // now create the truss and add it to the Domain |
| theTruss = new Truss2D(eleTag, iData[1], iData[2], *theMaterial, dData[0]); |
| if (theTruss == 0) { |
| opserr << "WARNING ran out of memory creating element with tag " << eleTag << endln; |
| delete theMaterial; |
| return 0; |
| } |
| return theTruss; |
| } |
| // typical constructor |
| Truss2D::Truss2D(int tag, |
| int Nd1, int Nd2, |
| UniaxialMaterial &theMat, |
| double a) |
| :Element(tag, 0), |
| externalNodes(2), |
| trans(1,4), L(0.0), A(a) |
| { |
| // get a copy of the material object for our own use |
| theMaterial = theMat.getCopy(); |
| if (theMaterial == 0) { |
| opserr << "FATAL Truss2D::Truss2D() - out of memory, could not get a copy of the Material\n"; |
| exit(-1); |
| } |
| // fill in the ID containing external node info with node id's |
| if (externalNodes.Size() != 2) { |
| opserr << "FATAL Truss2D::Truss2D() - out of memory, could not create an ID of size 2\n"; |
| exit(-1); |
| } |
| externalNodes(0) = Nd1; |
| externalNodes(1) = Nd2; |
| theNodes[0] = 0; |
| theNodes[1] = 0; |
| } |
| // constructor which should be invoked by an FE_ObjectBroker only |
| Truss2D::Truss2D() |
| :Element(0, 0), |
| theMaterial(0), |
| externalNodes(2), |
| trans(1,4), L(0.0), A(0.0) |
| { |
| theNodes[0] = 0; |
| theNodes[1] = 0; |
| } |
| // destructor - provided to clean up any memory |
| Truss2D::~Truss2D() |
| { |
| // clean up the memory associated with the element, this is |
| // memory the Truss2D objects allocates and memory allocated |
| // by other objects that the Truss2D object is responsible for |
| // cleaning up, i.e. the MaterialObject. |
| if (theMaterial != 0) |
| delete theMaterial; |
| } |
| int |
| Truss2D::getNumExternalNodes(void) const |
| { |
| return 2; |
| } |
| const ID & |
| Truss2D::getExternalNodes(void) |
| { |
| return externalNodes; |
| } |
| Node ** |
| Truss2D::getNodePtrs(void) |
| { |
| return theNodes; |
| } |
| int |
| Truss2D::getNumDOF(void) { |
| return 4; |
| } |
| // method: setDomain() |
| // to set a link to the enclosing Domain, ensure nodes exist in Domain |
| // and set pointers to these nodes, also determines the length and |
| // transformation Matrix. |
| void |
| Truss2D::setDomain(Domain *theDomain) |
| { |
| // check Domain is not null - invoked when object removed from a domain |
| if (theDomain == 0) { |
| return; |
| } |
| // first ensure nodes exist in Domain and set the node pointers |
| Node *end1Ptr, *end2Ptr; |
| int Nd1 = externalNodes(0); |
| int Nd2 = externalNodes(1); |
| end1Ptr = theDomain->getNode(Nd1); |
| end2Ptr = theDomain->getNode(Nd2); |
| if (end1Ptr == 0) { |
| opserr << "WARNING Truss2D::setDomain() - at truss " << this->getTag() << " node " << |
| Nd1 << " does not exist in domain\n"; |
| return; // don't go any further - otherwise segemntation fault |
| } |
| if (end2Ptr == 0) { |
| opserr << "WARNING Truss2D::setDomain() - at truss " << this->getTag() << " node " << |
| Nd2 << " does not exist in domain\n"; |
| return; // don't go any further - otherwise segemntation fault |
| } |
| theNodes[0] = end1Ptr; |
| theNodes[1] = end2Ptr; |
| // call the DomainComponent class method THIS IS VERY IMPORTANT |
| this->DomainComponent::setDomain(theDomain); |
| // ensure connected nodes have correct number of dof's |
| int dofNd1 = end1Ptr->getNumberDOF(); |
| int dofNd2 = end2Ptr->getNumberDOF(); |
| if ((dofNd1 != 2) || (dofNd2 != 2)) { |
| opserr << "Truss2D::setDomain(): 2 dof required at nodes\n"; |
| return; |
| } |
| // now determine the length & transformation matrix |
| const Vector &end1Crd = end1Ptr->getCrds(); |
| const Vector &end2Crd = end2Ptr->getCrds(); |
| double dx = end2Crd(0)-end1Crd(0); |
| double dy = end2Crd(1)-end1Crd(1); |
| L = sqrt(dx*dx + dy*dy); |
| if (L == 0.0) { |
| opserr << "WARNING Truss2D::setDomain() - Truss2D " << this->getTag() << |
| " has zero length\n"; |
| return; // don't go any further - otherwise divide by 0 error |
| } |
| double cs = dx/L; |
| double sn = dy/L; |
| trans(0,0) = -cs; |
| trans(0,1) = -sn; |
| trans(0,2) = cs; |
| trans(0,3) = sn; |
| } |
| int |
| Truss2D::commitState() |
| { |
| return theMaterial->commitState(); |
| } |
| int |
| Truss2D::revertToLastCommit() |
| { |
| return theMaterial->revertToLastCommit(); |
| } |
| int |
| Truss2D::revertToStart() |
| { |
| return theMaterial->revertToStart(); |
| } |
| int |
| Truss2D::update() |
| { |
| // determine the current strain given trial displacements at nodes |
| double strain = this->computeCurrentStrain(); |
| // set the strain in the materials |
| theMaterial->setTrialStrain(strain); |
| return 0; |
| } |
| const Matrix & |
| Truss2D::getTangentStiff(void) |
| { |
| if (L == 0.0) { // length = zero - problem in setDomain() warning message already printed |
| trussK.Zero(); |
| return trussK; |
| } |
| // get the current E from the material for the last updated strain |
| double E = theMaterial->getTangent(); |
| // form the tangent stiffness matrix |
| trussK = trans^trans; |
| trussK *= A*E/L; |
| // return the matrix |
| return trussK; |
| } |
| const Matrix & |
| Truss2D::getInitialStiff(void) |
| { |
| if (L == 0.0) { // length = zero - problem in setDomain() warning message already printed |
| trussK.Zero(); |
| return trussK; |
| } |
| // get the current E from the material for the last updated strain |
| double E = theMaterial->getInitialTangent(); |
| // form the tangent stiffness matrix |
| trussK = trans^trans; |
| trussK *= A*E/L; |
| // return the matrix |
| return trussK; |
| } |
| const Vector & |
| Truss2D::getResistingForce() |
| { |
| if (L == 0.0) { // if length == 0, problem in setDomain() |
| trussR.Zero(); |
| return trussR; |
| } |
| // want: R = Ku - Pext |
| // force = F * transformation |
| double force = A*theMaterial->getStress(); |
| for (int i=0; i<4; i++) |
| trussR(i) = trans(0,i)*force; |
| return trussR; |
| } |
| int |
| Truss2D::sendSelf(int commitTag, Channel &theChannel) |
| { |
| int res; |
| // note: we don't check for dataTag == 0 for Element |
| // objects as that is taken care of in a commit by the Domain |
| // object - don't want to have to do the check if sending data |
| int dataTag = this->getDbTag(); |
| // Truss2D packs it's data into a Vector and sends this to theChannel |
| // along with it's dbTag and the commitTag passed in the arguments |
| Vector data(5); |
| data(0) = this->getTag(); |
| data(1) = A; |
| data(2) = theMaterial->getClassTag(); |
| int matDbTag = theMaterial->getDbTag(); |
| // NOTE: we do have to ensure that the material has a database |
| // tag if we are sending to a database channel. |
| if (matDbTag == 0) { |
| matDbTag = theChannel.getDbTag(); |
| if (matDbTag != 0) |
| theMaterial->setDbTag(matDbTag); |
| } |
| data(3) = matDbTag; |
| res = theChannel.sendVector(dataTag, commitTag, data); |
| if (res < 0) { |
| opserr << "WARNING Truss2D::sendSelf() - failed to send Vector\n"; |
| return -1; |
| } |
| // Truss2D then sends the tags of it's two end nodes |
| res = theChannel.sendID(dataTag, commitTag, externalNodes); |
| if (res < 0) { |
| opserr << "WARNING Truss2D::sendSelf() - failed to send ID\n"; |
| return -2; |
| } |
| // finally Truss2D asks it's material object to send itself |
| res = theMaterial->sendSelf(commitTag, theChannel); |
| if (res < 0) { |
| opserr << "WARNING Truss2D::sendSelf() - failed to send the Material\n"; |
| return -3; |
| } |
| return 0; |
| } |
| int |
| Truss2D::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) |
| { |
| int res; |
| int dataTag = this->getDbTag(); |
| // Truss2D creates a Vector, receives the Vector and then sets the |
| // internal data with the data in the Vector |
| Vector data(5); |
| res = theChannel.recvVector(dataTag, commitTag, data); |
| if (res < 0) { |
| opserr << "WARNING Truss2D::recvSelf() - failed to receive Vector\n"; |
| return -1; |
| } |
| this->setTag((int)data(0)); |
| A = data(1); |
| // Truss2D now receives the tags of it's two external nodes |
| res = theChannel.recvID(dataTag, commitTag, externalNodes); |
| if (res < 0) { |
| opserr << "WARNING Truss2D::recvSelf() - failed to receive ID\n"; |
| return -2; |
| } |
| // we create a material object of the correct type, |
| // sets its database tag and asks this new object to recveive itself. |
| int matClass = data(2); |
| int matDb = data(3); |
| theMaterial = theBroker.getNewUniaxialMaterial(matClass); |
| if (theMaterial == 0) { |
| opserr << "WARNING Truss2D::recvSelf() - failed to create a Material\n"; |
| return -3; |
| } |
| // we set the dbTag before we receive the material - this is important |
| theMaterial->setDbTag(matDb); |
| res = theMaterial->recvSelf(commitTag, theChannel, theBroker); |
| if (res < 0) { |
| opserr << "WARNING Truss2D::recvSelf() - failed to receive the Material\n"; |
| return -3; |
| } |
| return 0; |
| } |
| void |
| Truss2D::Print(OPS_Stream &s, int flag) |
| { |
| s << "Element: " << this->getTag(); |
| s << " type: Truss2D iNode: " << externalNodes(0); |
| s << " jNode: " << externalNodes(1); |
| s << " Area: " << A; |
| s << " \t Material: " << *theMaterial; |
| } |
| Response * |
| Truss2D::setResponse(const char **argv, int argc, OPS_Stream &output) |
| { |
| Response *theResponse = 0; |
| output.tag("ElementOutput"); |
| output.attr("eleType",this->getClassType()); |
| output.attr("eleTag",this->getTag()); |
| int numNodes = this->getNumExternalNodes(); |
| const ID &nodes = this->getExternalNodes(); |
| static char nodeData[32]; |
| for (int i=0; i<numNodes; i++) { |
| sprintf(nodeData,"node%d",i+1); |
| output.attr(nodeData,nodes(i)); |
| } |
| if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0 || |
| strcmp(argv[0],"globalForce") == 0 || strcmp(argv[0],"globalForces") == 0) { |
| const Vector &force = this->getResistingForce(); |
| int size = force.Size(); |
| for (int i=0; i<size; i++) { |
| sprintf(nodeData,"P%d",i+1); |
| output.tag("ResponseType",nodeData); |
| } |
| theResponse = new ElementResponse(this, 1, this->getResistingForce()); |
| } |
| else if (strcmp(argv[0],"dampingForce") == 0 || strcmp(argv[0],"dampingForces") == 0) { |
| const Vector &force = this->getResistingForce(); |
| int size = force.Size(); |
| for (int i=0; i<size; i++) { |
| sprintf(nodeData,"P%d",i+1); |
| output.tag("ResponseType",nodeData); |
| } |
| theResponse = new ElementResponse(this, 2, this->getResistingForce()); |
| } else if (strcmp(argv[0],"axialForce") ==0) |
| return new ElementResponse(this, 3, 0.0); |
| output.endTag(); |
| return theResponse; |
| } |
| int |
| Truss2D::getResponse(int responseID, Information &eleInfo) |
| { |
| double strain; |
| switch (responseID) { |
| case -1: |
| return -1; |
| case 1: // global forces |
| return eleInfo.setVector(this->getResistingForce()); |
| case 2: |
| return eleInfo.setVector(this->getRayleighDampingForces()); |
| case 3: |
| theMaterial->setTrialStrain(strain); |
| return eleInfo.setDouble(A * theMaterial->getStress()); |
| default: |
| return 0; |
| } |
| } |
| double |
| Truss2D::computeCurrentStrain(void) const |
| { |
| // NOTE this method will never be called with L == 0.0 |
| // determine the strain |
| const Vector &disp1 = theNodes[0]->getTrialDisp(); |
| const Vector &disp2 = theNodes[1]->getTrialDisp(); |
| double dLength = 0.0; |
| for (int i=0; i<2; i++){ |
| dLength -= (disp2(i)-disp1(i)) * trans(0,i); |
| } |
| double strain = dLength/L; |
| return strain; |
| } |
| /trunk/DEVELOPER/element/cpp/Truss2D.h |
|---|
| New file |
| 0,0 → 1,106 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision: 1.2 $ |
| // $Date: 2008/12/10 00:05:21 $ |
| // $Source: /usr/local/cvs/OpenSees/PACKAGES/NewElement/cpp/Truss2D.h,v $ |
| #ifndef Truss2D_h |
| #define Truss2D_h |
| // Written: fmk |
| // |
| // Description: This file contains the interface for the Truss2D class. |
| // It defines the class interface and the class attributes. Truss2D |
| // provides the abstraction of a simple truss element for 2-d problems. |
| // The stress-strain relationship for the truss being performed in the |
| // associated UniaxialMaterial object. |
| // |
| // What: "@(#) Truss2D.h, revA" |
| #include <Element.h> |
| #include <Matrix.h> |
| #include <Vector.h> |
| class UniaxialMaterial; |
| class Truss2D : public Element |
| { |
| public: |
| // constructors |
| Truss2D(int tag, |
| int Nd1, int Nd2, |
| UniaxialMaterial &theMaterial, |
| double A); |
| Truss2D(); |
| // destructor |
| ~Truss2D(); |
| // public methods to obtain inforrmation about dof & connectivity |
| int getNumExternalNodes(void) const; |
| const ID &getExternalNodes(void); |
| Node **getNodePtrs(void); |
| int getNumDOF(void); |
| void setDomain(Domain *theDomain); |
| // public methods to set the state of the element |
| int commitState(void); |
| int revertToLastCommit(void); |
| int revertToStart(void); |
| int update(void); |
| // public methods to obtain stiffness |
| const Matrix &getTangentStiff(void); |
| const Matrix &getInitialStiff(void); |
| // public method to obtain resisting force |
| const Vector &getResistingForce(void); |
| // public methods for output |
| int sendSelf(int commitTag, Channel &theChannel); |
| int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker); |
| void Print(OPS_Stream &s, int flag =0); |
| Response *setResponse(const char **argv, int argc, OPS_Stream &s); |
| int getResponse(int responseID, Information &eleInformation); |
| protected: |
| private: |
| // private member functions - only available to objects of the class |
| double computeCurrentStrain(void) const; |
| // private attributes - a copy for each object of the class |
| UniaxialMaterial *theMaterial; // pointer to a material |
| ID externalNodes; // contains the id's of end nodes |
| Matrix trans; // hold the transformation matrix |
| double L; // length of truss (undeformed configuration) - set in setDomain() |
| double A; // area of truss |
| Node *theNodes[2]; // node pointers |
| // static data - single copy for all objects of the class |
| static Matrix trussK; // class wide matrix for returning stiffness |
| static Vector trussR; // class wide vector for returning residual |
| }; |
| #endif |
| /trunk/DEVELOPER/element/cpp/Makefile |
|---|
| New file |
| 0,0 → 1,41 |
| include ../../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| all: $(LIBRARY_CPP) $(LIBRARY_C) $(LIBRARY_F) |
| LIBRARY_CPP = Truss2D.so |
| OBJS_CPP = Truss2D.o |
| $(LIBRARY_CPP): $(OBJS_CPP) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY_CPP) $(OBJS_CPP) -lc |
| Truss2D.o: Truss2D.cpp |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall Truss2D.cpp |
| else |
| LIBRARY_CPP = Truss2D.dylib |
| all: $(OBJS) |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace Truss2D.cpp $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY_CPP) |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/element/cpp/example1.tcl |
|---|
| New file |
| 0,0 → 1,58 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.2 $ |
| # $Date: 2008/12/09 22:43:40 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewElement/cpp/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial Elastic 1 3000 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element Truss2D 1 1 4 10.0 1 |
| element Truss2D 2 2 4 5.0 1 |
| element Truss2D 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system BandSPD |
| constraints Plain |
| integrator LoadControl 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| exit |
| /trunk/DEVELOPER/element/fortran/elementTypes.f |
|---|
| New file |
| 0,0 → 1,62 |
| module elementTypes |
| use iso_c_binding |
| implicit none |
| private |
| public modelState |
| type, bind(c) :: modelState |
| real (C_DOUBLE) time |
| real (C_DOUBLE) dt |
| end type modelState |
| c type modelState |
| c real *8::time; |
| c real *8::dt; |
| c end type modelState |
| public matObject |
| type, bind(C) :: matObject |
| integer(C_INT) ::tag; |
| integer(C_INT) ::nParam; |
| integer(C_INT) ::nState |
| type (c_ptr) :: theParam; |
| type (c_ptr) :: cState; |
| type (c_ptr) :: tState; |
| type(c_funptr) ::functPtr; |
| type (c_ptr) :: matObjPtr; |
| end type matObject; |
| c |
| c real (C_DOUBLE) :: theParam; |
| c real (C_DOUBLE) :: cState; |
| c real (C_DOUBLE) :: tState; |
| c public matObject |
| c type :: matObject |
| c integer::tag |
| c integer::nParam |
| c integer::nState |
| c real *8, pointer ::theParam(:) |
| c real *8, pointer ::cState(:) |
| c real *8, pointer ::tState(:) |
| c type(c_funptr)::functPtr |
| c end type matObject |
| public eleObject |
| type, bind(C) :: eleObject |
| integer(C_INT) :: tag; |
| integer(C_INT) :: nNode; |
| integer(C_INT) :: nDOF; |
| integer(C_INT) :: nParam; |
| integer(C_INT) :: nState; |
| integer(C_INT) :: nMat; |
| type (c_ptr) :: node; |
| type (c_ptr) :: param; |
| type (c_ptr) :: cState; |
| type (c_ptr) :: tState; |
| type (c_ptr) :: mats; |
| type(c_funptr) ::functPtr |
| end type eleObject; |
| end module |
| /trunk/DEVELOPER/element/fortran/elementAPI.f |
|---|
| New file |
| 0,0 → 1,168 |
| module elementAPI |
| use iso_c_binding |
| use elementTypes |
| implicit none |
| integer ::ISW_COMMIT, ISW_INIT, ISW_REVERT_TO_START |
| integer ::ISW_REVERT, ISW_FORM_TANG_AND_RESID, ISW_FORM_MASS |
| integer ::OPS_UNIAXIAL_MATERIAL_TYPE |
| real *8 ::DBL_EPSILON |
| PARAMETER (ISW_INIT = 0) |
| PARAMETER (ISW_COMMIT = 1) |
| PARAMETER (ISW_REVERT = 2) |
| PARAMETER (ISW_FORM_TANG_AND_RESID = 3) |
| PARAMETER (ISW_FORM_MASS = 4) |
| PARAMETER (ISW_REVERT_TO_START = 5) |
| PARAMETER (DBL_EPSILON = 1.0e-21) |
| PARAMETER (OPS_UNIAXIAL_MATERIAL_TYPE = 1) |
| PUBLIC OPS_GetIntInput |
| interface |
| function OPS_GetIntInput(numData, iData) |
| integer :: OPS_GetIntInput |
| integer :: numData |
| integer :: iData(*) |
| end function OPS_GetIntInput |
| end interface |
| PUBLIC OPS_GetDoubleInput |
| interface |
| function OPS_GetDoubleInput(numData, dData) |
| integer :: OPS_GetDoubleInput |
| integer :: numData |
| real *8 :: dData(*) |
| end function OPS_GetDoubleInput |
| end interface |
| PUBLIC OPS_AllocateMaterial |
| interface |
| function OPS_AllocateMaterial(mat) |
| use elementtypes |
| integer :: OPS_AllocateMaterial |
| type(matObject) :: mat |
| end function OPS_AllocateMaterial |
| end interface |
| PUBLIC OPS_GetMaterial |
| interface |
| function OPS_GetMaterial(matTag, matType) |
| use iso_c_binding |
| use elementtypes |
| type(c_ptr) :: OPS_GetMaterial |
| integer :: matTag |
| integer :: matType |
| end function OPS_GetMaterial |
| end interface |
| c integer :: matTags(*) |
| c type(c_ptr) :: matTags |
| PUBLIC OPS_AllocateElement |
| interface |
| function OPS_AllocateElement(ele, matTags, matType) |
| use elementtypes |
| use iso_c_binding |
| integer :: OPS_AllocateElement |
| type(eleObject) :: ele |
| integer :: matTags(*) |
| integer :: matType |
| end function OPS_AllocateElement |
| end interface |
| PUBLIC OPS_GetNodeCrd |
| interface |
| function OPS_GetNodeCrd(nodeTag, numData, dData) |
| integer :: OPS_GetNodeCrd |
| integer :: nodeTag |
| integer :: numData |
| real *8 :: dData(*) |
| end function OPS_GetNodeCrd |
| end interface |
| PUBLIC OPS_GetNodeDisp |
| interface |
| function OPS_GetNodeDisp(nodeTag, numData, dData) |
| integer :: OPS_GetNodeDisp |
| integer :: nodeTag |
| integer :: numData |
| real *8 :: dData(*) |
| end function OPS_GetNodeDisp |
| end interface |
| PUBLIC OPS_GetNodeIncrDisp |
| interface |
| function OPS_GetNodeIncrDisp(nodeTag, numData, dData) |
| integer :: OPS_GetNodeDisp |
| integer :: nodeTag |
| integer :: numData |
| real *8 :: dData(*) |
| end function OPS_GetNodeIncrDisp |
| end interface |
| PUBLIC OPS_GetNodeIncrDeltaDisp |
| interface |
| function OPS_GetNodeIncrDeltaDisp(nodeTag, numData, dData) |
| integer :: OPS_GetNodeDisp |
| integer :: nodeTag |
| integer :: numData |
| real *8 :: dData(*) |
| end function OPS_GetNodeIncrDeltaDisp |
| end interface |
| PUBLIC OPS_InvokeMaterial |
| interface |
| function OPS_InvokeMaterial(e, mt, md, strn, strs, tang, isw) |
| use elementtypes |
| integer :: OPS_InvokeMaterial |
| type(eleObject) :: e |
| integer :: mt |
| type(modelState) :: md |
| real *8 :: strn(*) |
| real *8 :: strs(*) |
| real *8 :: tang(*) |
| integer :: isw |
| end function OPS_InvokeMaterial |
| end interface |
| PUBLIC OPS_InvokeMaterialDirectly |
| interface |
| function OPS_InvokeMaterialDirectly(mat, md, strn, strs, tang, i) |
| use elementtypes |
| use iso_c_binding |
| integer :: OPS_InvokeMaterialDirectly |
| type(c_ptr) :: mat |
| type(modelState) :: md |
| real *8 :: strn(*) |
| real *8 :: strs(*) |
| real *8 :: tang(*) |
| integer :: i |
| end function OPS_InvokeMaterialDirectly |
| end interface |
| PUBLIC OPS_InvokeMaterialDirectly2 |
| interface |
| function OPS_InvokeMaterialDirectly2(mat, md, strn, strs, tang, i) |
| use elementtypes |
| use iso_c_binding |
| integer :: OPS_InvokeMaterialDirectly2 |
| type(c_ptr), value :: mat |
| type(modelState) :: md |
| real *8 :: strn(*) |
| real *8 :: strs(*) |
| real *8 :: tang(*) |
| integer :: i |
| end function OPS_InvokeMaterialDirectly2 |
| end interface |
| PUBLIC OPS_Error |
| interface |
| function OPS_Error(msg) |
| integer :: OPS_Error |
| character, dimension(*) :: msg |
| end function OPS_Error |
| end interface |
| end module elementAPI |
| /trunk/DEVELOPER/element/fortran/trussF.f |
|---|
| New file |
| 0,0 → 1,198 |
| 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(4, *) |
| double precision resid(1) |
| 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 |
| 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 |
| 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 = 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 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; |
| do 20 i =1,4 |
| resid(i) = force * tran(i); |
| do 30 j = 1,4 |
| tang(i,j) = k * tran(i)*tran(j); |
| 30 continue |
| 20 continue |
| END IF |
| END IF |
| c return error code |
| error = 0 |
| END SUBROUTINE trussf |
| SUBROUTINE localinit() |
| !DEC$ IF DEFINED (_DLL) |
| !DEC$ ATTRIBUTES DLLEXPORT :: LOCALINIT |
| !DEC$ END IF |
| use elementAPI |
| implicit none |
| integer::error; |
| character *60:: msg; |
| msg = 'trussf - Written by fmk UC Berkeley Copyright 2008' |
| error = OPS_Error(msg); |
| END SUBROUTINE localInit |
| /trunk/DEVELOPER/element/fortran/Makefile |
|---|
| New file |
| 0,0 → 1,46 |
| include ../../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| all: $(LIBRARY_F) |
| LIBRARY_F = trussf.so |
| OBJS_F = elementtypes.mod elementapi.mod trussF.o |
| $(LIBRARY_F): $(OBJS_F) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY) trussF.o -lc |
| else |
| LIBRARY_F = trussf.dylib |
| OTHER_OBJS = dynamicLoadFunctions.o |
| OBJS = elementtypes.mod elementapi.mod |
| all: $(OBJS) |
| $(FC) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace trussf.f $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY_F) |
| elementapi.mod: elementAPI.f elementtypes.mod |
| $(FC) -c elementAPI.f |
| elementtypes.mod: elementTypes.f |
| $(FC) -c elementTypes.f |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/element/fortran/example1.tcl |
|---|
| New file |
| 0,0 → 1,57 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.2 $ |
| # $Date: 2009/01/13 07:31:06 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewElement/fortran/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.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 4 10.0 1 |
| element trussf 2 2 4 5.0 1 |
| element trussf 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system BandSPD |
| constraints Plain |
| integrator LoadControl 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| /trunk/DEVELOPER/element/Makefile |
|---|
| New file |
| 0,0 → 1,24 |
| include ../Makefile.def |
| all: |
| echo @pwd; |
| @$(CD) cpp; $(MAKE); |
| @$(CD) c; $(MAKE); |
| @$(CD) fortran; $(MAKE); |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(CD) ./cpp; $(MAKE) wipe; |
| @$(CD) ./c; $(MAKE) wipe; |
| @$(CD) ./fortran; $(MAKE) wipe; |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| wipe: spotless |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/Makefile.def |
|---|
| New file |
| 0,0 → 1,30 |
| OS_FLAG = -D_MACOSX |
| #OS_FLAG = -D_LINUX |
| CC++ = g++ |
| CC = gcc |
| FC = /usr/local/bin/gfortran |
| LINKER = $(CC++) |
| LINKFLAGS = -Wl |
| CD = cd |
| C++FLAGS = -O3 |
| CFLAGS = -O2 |
| FFLAGS = -O |
| INCLUDES = -I../core -I../../core |
| .cpp.o: |
| @$(ECHO) Making $@ from $< |
| $(CC++) $(C++FLAGS) $(INCLUDES) -c $< -o $@ |
| .C.o: |
| @$(ECHO) Making $@ from $< |
| $(CC++) $(C++FLAGS) $(INCLUDES) -c $< -o $@ |
| .c.o: |
| @$(ECHO) Making $@ from $< |
| $(CC) $(CFLAGS) -c $< -o $@ |
| .f.o: |
| @$(ECHO) Making $@ from $< |
| $(FC) $(FFLAGS) -c $< -o $@ |
| /trunk/DEVELOPER/Makefile |
|---|
| New file |
| 0,0 → 1,25 |
| include ./Makefile.def |
| all: |
| @cd ./element; $(MAKE); |
| @cd ./material; $(MAKE); |
| @cd ./recorder; $(MAKE); |
| @cd ./system; $(MAKE); |
| @cd ./integrator; $(MAKE); |
| tidy: |
| @rm -f Makefile.bak *~ #*# core |
| clean: tidy |
| @rm -f *.o |
| spotless: clean |
| wipe: spotless |
| @cd ./element; $(MAKE) wipe; |
| @cd ./material; $(MAKE) wipe; |
| @cd ./recorder; $(MAKE) wipe; |
| @cd ./system; $(MAKE) wipe; |
| @cd ./integrator; $(MAKE) wipe; |
| /trunk/DEVELOPER/WindowsInstructions.txt |
|---|
| New file |
| 0,0 → 1,45 |
| 1. Open ViualStudio |
| 2. File -> New Project |
| 3. Give it name of New class and location of current class files (ElasticPPCPP and C:\???\OpenSeesDeveloper\material |
| 4. Select OK. |
| 5. New windows pops up, Select Application Settins |
| 6. Select DLL AND select Empty Project |
| 7. Select Finish |
| a new project pops up in the workspace with class name (ElasticPPCPP) |
| 8. Right click on source files, select add Existing Item. |
| 9. browse to the class file (should be up 1 directory, select it (elasticPPcPP.cpp) |
| 10.Right click on heared files, select add existing item. |
| 11. browse to header file, select it (elasticOOCPP.h) |
| 12. right click on Project (ElasticPPCPP), select build. |
| It fails in compilation, as cannot find file elementAPI.h |
| 13. right click on Project (ElasticPPCPP), select Properties. |
| 14. select C/C++ folder icon |
| 15. In Configuration pull down menu, select all Configurations |
| 16. Click 3 dots to right of Additional Include Directories |
| 17. In window that pops up, select folder |
| 18. add to line ..\..\core (this directory contains the elementAPI.h file). |
| 19. select ok. |
| 20 right click on project (ElasticPPCPP), select build |
| It fails in Linking, a lot of unresolved external symbols. |
| 21. Right click on source files, select add Existing Item. |
| 22. browse to the core directory (should be 1 directory up) |
| 23. select all the .cpp files here |
| 24. right click on Project (ElasticPPCPP), select build. |
| IT SHOULD WORK! .. IF IT FAILS GET MY ATTENTION. |
| Now test it, |
| 25. using windows finder, copy the .dll (C:\??\OpenSeesDeveloper\material\ElasticPPCPP\Debug) into the class directory (C:\???\material\) |
| 26. run OpenSees (2.2.0) (C:\??\OpenSeesDeveloper\bin) |
| 27. cd into directory containing example (cd ..\material) |
| 28. source an example script containing your new command, i.e. type: "source example1.tcl" |
| IT SHOULD RUN THE SCRIPT! .. IF IT FAILS GET MY ATTENTION. |
| /trunk/DEVELOPER/integrator/Trapezoidal.cpp |
|---|
| New file |
| 0,0 → 1,379 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision: 1.1 $ |
| // $Date: 2009/03/20 22:37:09 $ |
| // $Source: /usr/local/cvs/OpenSees/PACKAGES/NewIntegrator/Trapezoidal.cpp,v $ |
| // Written : fmk |
| // Created : 02/09 |
| // |
| // Description: This file contains the implementation of the Trapezoidal class. |
| // ref: K.J.Bathe, "Conserving Energy and Momentum in Nonlinear Dynamics: A Simple |
| // Implicit Time Integration Scheme", Computers ans Structures 85(2007),437-445 |
| // |
| // NOTE: In the implementation we use deltaT as opposed to deltaT/2 in the paper. |
| // |
| // What: "@(#) Trapezoidal.C, revA" |
| #include <elementAPI.h> |
| #include "Trapezoidal.h" |
| #include <FE_Element.h> |
| #include <LinearSOE.h> |
| #include <AnalysisModel.h> |
| #include <Vector.h> |
| #include <DOF_Group.h> |
| #include <DOF_GrpIter.h> |
| #include <AnalysisModel.h> |
| #include <Channel.h> |
| #include <FEM_ObjectBroker.h> |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export extern "C" _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| // init function called once on first loading of package |
| OPS_Export void |
| localInit() |
| { |
| OPS_Error("Trapezoidal - Written by fmk UC Berkeley Copyright 2009\n", 1); |
| } |
| OPS_Export void * |
| OPS_Trapezoidal() |
| { |
| // |
| // parse args |
| // |
| // none to parse for this soe |
| // |
| // create the SOE |
| // |
| TransientIntegrator *theIntegrator = new Trapezoidal(); |
| // return it |
| return theIntegrator; |
| } |
| Trapezoidal::Trapezoidal() |
| : TransientIntegrator(0), |
| c1(0.0), c2(0.0), c3(0.0), |
| Ut(0), Utdot(0), Utdotdot(0), |
| U(0), Udot(0), Udotdot(0) |
| { |
| } |
| Trapezoidal::~Trapezoidal() |
| { |
| if (Ut != 0) |
| delete Ut; |
| if (Utdot != 0) |
| delete Utdot; |
| if (Utdotdot != 0) |
| delete Utdotdot; |
| if (U != 0) |
| delete U; |
| if (Udot != 0) |
| delete Udot; |
| if (Udotdot != 0) |
| delete Udotdot; |
| } |
| int Trapezoidal::newStep(double deltaT) |
| { |
| // check the vectors have been created |
| if (U == 0) { |
| opserr << "Trapezoidal::newStep() - domainChange() failed or hasn't been called\n"; |
| return -3; |
| } |
| // get a pointer to the AnalysisModel |
| AnalysisModel *theModel = this->getAnalysisModel(); |
| // set response at t to be that at t+deltaT of previous step |
| (*Ut) = *U; |
| (*Utdot) = *Udot; |
| (*Utdotdot) = *Udotdot; |
| c1 = 1.0; |
| c2 = 2.0/deltaT; |
| c3 = 4.0/(deltaT*deltaT); |
| (*Udot) *= -1.0; |
| double a3 = -4.0/deltaT; |
| double a4 = -1; |
| Udotdot->addVector(a4, *Utdot, a3); |
| // set the trial response quantities |
| theModel->setVel(*Udot); |
| theModel->setAccel(*Udotdot); |
| // increment the time to t+deltaT and apply the load |
| double time = theModel->getCurrentDomainTime(); |
| time += deltaT; |
| if (theModel->updateDomain(time, deltaT) < 0) { |
| opserr << "Trapezoidal::newStep() - failed to update the domain\n"; |
| return -4; |
| } |
| return 0; |
| } |
| int Trapezoidal::revertToLastStep() |
| { |
| // set response at t+deltaT to be that at t .. for next newStep |
| if (U != 0) { |
| (*U) = *Ut; |
| (*Udot) = *Utdot; |
| (*Udotdot) = *Utdotdot; |
| } |
| return 0; |
| } |
| int Trapezoidal::formEleTangent(FE_Element *theEle) |
| { |
| theEle->zeroTangent(); |
| if (statusFlag == CURRENT_TANGENT) { |
| theEle->addKtToTang(c1); |
| theEle->addCtoTang(c2); |
| theEle->addMtoTang(c3); |
| } else if (statusFlag == INITIAL_TANGENT) { |
| theEle->addKiToTang(c1); |
| theEle->addCtoTang(c2); |
| theEle->addMtoTang(c3); |
| } |
| return 0; |
| } |
| int Trapezoidal::formNodTangent(DOF_Group *theDof) |
| { |
| theDof->zeroTangent(); |
| theDof->addCtoTang(c2); |
| theDof->addMtoTang(c3); |
| return 0; |
| } |
| int Trapezoidal::domainChanged() |
| { |
| AnalysisModel *myModel = this->getAnalysisModel(); |
| LinearSOE *theLinSOE = this->getLinearSOE(); |
| const Vector &x = theLinSOE->getX(); |
| int size = x.Size(); |
| // create the new Vector objects |
| if (Ut == 0 || Ut->Size() != size) { |
| if (Ut != 0) |
| delete Ut; |
| if (Utdot != 0) |
| delete Utdot; |
| if (Utdotdot != 0) |
| delete Utdotdot; |
| if (U != 0) |
| delete U; |
| if (Udot != 0) |
| delete Udot; |
| if (Udotdot != 0) |
| delete Udotdot; |
| // create the new |
| Ut = new Vector(size); |
| Utdot = new Vector(size); |
| Utdotdot = new Vector(size); |
| U = new Vector(size); |
| Udot = new Vector(size); |
| Udotdot = new Vector(size); |
| // check we obtained the new |
| if (Ut == 0 || Ut->Size() != size || |
| Utdot == 0 || Utdot->Size() != size || |
| Utdotdot == 0 || Utdotdot->Size() != size || |
| U == 0 || U->Size() != size || |
| Udot == 0 || Udot->Size() != size || |
| Udotdot == 0 || Udotdot->Size() != size) { |
| // delete any we obtained |
| if (Ut != 0) |
| delete Ut; |
| if (Utdot != 0) |
| delete Utdot; |
| if (Utdotdot != 0) |
| delete Utdotdot; |
| if (U != 0) |
| delete U; |
| if (Udot != 0) |
| delete Udot; |
| if (Udotdot != 0) |
| delete Udotdot; |
| Ut = 0; Utdot = 0; Utdotdot = 0; |
| U = 0; Udot = 0; Udotdot = 0; |
| return -1; |
| } |
| } |
| // now go through and populate U, Udot and Udotdot by iterating through |
| // the DOF_Groups and getting the last committed velocity and accel |
| DOF_GrpIter &theDOFs = myModel->getDOFs(); |
| DOF_Group *dofPtr; |
| while ((dofPtr = theDOFs()) != 0) { |
| const ID &id = dofPtr->getID(); |
| int idSize = id.Size(); |
| int i; |
| const Vector &disp = dofPtr->getCommittedDisp(); |
| for (i=0; i < idSize; i++) { |
| int loc = id(i); |
| if (loc >= 0) { |
| (*U)(loc) = disp(i); |
| } |
| } |
| const Vector &vel = dofPtr->getCommittedVel(); |
| for (i=0; i < idSize; i++) { |
| int loc = id(i); |
| if (loc >= 0) { |
| (*Udot)(loc) = vel(i); |
| } |
| } |
| const Vector &accel = dofPtr->getCommittedAccel(); |
| for (i=0; i < idSize; i++) { |
| int loc = id(i); |
| if (loc >= 0) { |
| (*Udotdot)(loc) = accel(i); |
| } |
| } |
| } |
| return 0; |
| } |
| int Trapezoidal::update(const Vector &deltaU) |
| { |
| AnalysisModel *theModel = this->getAnalysisModel(); |
| if (theModel == 0) { |
| opserr << "WARNING Trapezoidal::update() - no AnalysisModel set\n"; |
| return -1; |
| } |
| // check domainChanged() has been called, i.e. Ut will not be zero |
| if (Ut == 0) { |
| opserr << "WARNING Trapezoidal::update() - domainChange() failed or not called\n"; |
| return -2; |
| } |
| // check deltaU is of correct size |
| if (deltaU.Size() != U->Size()) { |
| opserr << "WARNING Trapezoidal::update() - Vectors of incompatible size "; |
| opserr << " expecting " << U->Size() << " obtained " << deltaU.Size() << endln; |
| return -3; |
| } |
| (*U) += deltaU; |
| Udot->addVector(1.0, deltaU, c2); |
| Udotdot->addVector(1.0, deltaU, c3); |
| // update the response at the DOFs |
| theModel->setResponse(*U,*Udot,*Udotdot); |
| if (theModel->updateDomain() < 0) { |
| opserr << "Trapezoidal::update() - failed to update the domain\n"; |
| return -4; |
| } |
| return 0; |
| } |
| int Trapezoidal::sendSelf(int cTag, Channel &theChannel) |
| { |
| return 0; |
| } |
| int Trapezoidal::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker) |
| { |
| return 0; |
| } |
| void Trapezoidal::Print(OPS_Stream &s, int flag) |
| { |
| AnalysisModel *theModel = this->getAnalysisModel(); |
| if (theModel != 0) { |
| double currentTime = theModel->getCurrentDomainTime(); |
| s << "\t Trapezoidal - currentTime: " << currentTime; |
| } else |
| s << "\t Trapezoidal - no associated AnalysisModel\n"; |
| } |
| // AddingSensitivity:BEGIN ////////////////////////////// |
| int Trapezoidal::revertToStart() |
| { |
| if (Ut != 0) |
| Ut->Zero(); |
| if (Utdot != 0) |
| Utdot->Zero(); |
| if (Utdotdot != 0) |
| Utdotdot->Zero(); |
| if (U != 0) |
| U->Zero(); |
| if (Udot != 0) |
| Udot->Zero(); |
| if (Udotdot != 0) |
| Udotdot->Zero(); |
| return 0; |
| } |
| // AddingSensitivity:END //////////////////////////////// |
| /trunk/DEVELOPER/integrator/Trapezoidal.h |
|---|
| New file |
| 0,0 → 1,78 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision: 1.1 $ |
| // $Date: 2009/03/20 22:37:09 $ |
| // $Source: /usr/local/cvs/OpenSees/PACKAGES/NewIntegrator/Trapezoidal.h,v $ |
| #ifndef Trapezoidal_h |
| #define Trapezoidal_h |
| // Written : fmk |
| // Created : 02/09 |
| // |
| // Description: This file contains the class definition for Trapezoidal. |
| // Trapezoidal is a transient integrator which uses trapezoidal rule for |
| // Ut+deltaT and Vt+deltaT (U,V,A the displ, vel and accel) |
| // Ut+deltaT = Ut + [deltaT/2.0] * (Vt + Vt+deltaT) |
| // Vt+deltaT = Vt + [deltaT/2.0] * (At + At+deltaT) |
| // and solves transient equation at time t+deltaT. |
| // |
| // It is the numerically the same as Newmark with gamma=0.5, beta=0.25 (average acel method) |
| // What: "@(#) Trapezoidal.h, revA" |
| #include <TransientIntegrator.h> |
| class DOF_Group; |
| class FE_Element; |
| class Vector; |
| class Trapezoidal : public TransientIntegrator |
| { |
| public: |
| Trapezoidal(); |
| ~Trapezoidal(); |
| // methods which define what the FE_Element and DOF_Groups add |
| // to the system of equation object. |
| int formEleTangent(FE_Element *theEle); |
| int formNodTangent(DOF_Group *theDof); |
| int domainChanged(void); |
| int newStep(double deltaT); |
| int revertToLastStep(void); |
| int update(const Vector &deltaU); |
| int sendSelf(int commitTag, Channel &theChannel); |
| int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker); |
| void Print(OPS_Stream &s, int flag = 0); |
| int revertToStart(); |
| protected: |
| private: |
| double c1, c2, c3; // some constants we need to keep |
| Vector *Ut, *Utdot, *Utdotdot; // response quantities at time t |
| Vector *U, *Udot, *Udotdot; // response quantities at time t+deltaT |
| }; |
| #endif |
| /trunk/DEVELOPER/integrator/Makefile |
|---|
| New file |
| 0,0 → 1,40 |
| include ../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| LIBRARY = Trapezoidal.so |
| all: $(LIBRARY) |
| OBJS = Trapezoidal.o |
| $(LIBRARY): $(OBJS) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY) $(OBJS) -lc |
| TrussCPP.o: TrussCPP.cpp |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall TrussCPP.cpp |
| else |
| LIBRARY = Trapezoidal.dylib |
| all: |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace Trapezoidal.cpp $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY) |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/integrator/example1.tcl |
|---|
| New file |
| 0,0 → 1,158 |
| # Example: 1 element truss model subjcted to sine wave input motion at base |
| # |
| # model Uf,Af --> |______/\/\/\/\/\_____** M -->U,V,A |
| # | K ** |
| # |
| # subject to motion: Uf = USinwt |
| # |
| # | M 0 | A + | K -K|U = 0 |
| # | 0 M2| Af | -K K|Uf 0 |
| # |
| # rewriting eqn 1 of equtaions: |
| # |
| # MA + KU = KU Sinwt subject to: U(0) = 0, V(0) = 0; |
| # |
| # Exact Solution: |
| # |
| # U(t) = ____U_____ Sin(wt) - _U_ (w/wn)_ Sin(wnt) |
| # (1-(w/wn)^2) (1-(w/wn)^2) |
| # |
| # |
| # Uniform Excitation Case: Af=-Uw^2 Sinwt Vf(0) = -Uw |
| # |
| # U(t) = Uf(t) + Ur(t) |
| # |
| # MAr + KUr = -MAf = MUw^2 Sinwt, Ur(0) = 0; Vr(0) = -Uw |
| # |
| # Exact Solution for Uniform Case: |
| # |
| # Ur(t) = ___MUw^2____ Sin(wt) - Uw Sin(wnt) - _MUw^2 (w/wn)_ Sin(wnt) |
| # K(1-(w/wn)^2) wn K(1-(w/wn)^2) |
| # |
| # = _____Uw^2______ Sin(wt) - Uw Sin(wnt) - __Uw^2_ (w/wn)__ Sin(wnt) |
| # wn^2(1-(w/wn)^2) wn wn^2 (1-(w/wn)^2) |
| # |
| # = _____Uw^2______ Sin(wt) - ____ U w ____ Sin(wnt) |
| # wn^2(1-(w/wn)^2) wn(1-(w/wn)^2) |
| # |
| # U(t) = Ur(t) + Uf(t) |
| # |
| # U(t) = ____U_____ Sin(wt) - _U_ (w/wn)_ Sin(wnt) |
| # (1-(w/wn)^2) (1-(w/wn)^2) |
| # |
| proc doModel {K periodN U period tFinish inputT} { |
| # some constants |
| set g 386.4 |
| set PI 3.14159265 |
| # some clculated parameters for basic model |
| set omega [expr 2*$PI/$period] |
| set omegaN [expr 2*$PI/$periodN] |
| set M [expr $K/($omegaN*$omegaN)] |
| set factor $U |
| puts "$period $periodN $omega $omegaN" |
| # |
| # Create Model |
| # |
| model BasicBuilder -ndm 1 -ndf 1; # Define the model builder, ndm=#dimension, ndf=#dofs |
| wipe |
| node 1 0 # node#, X, Y |
| node 2 1.0 -mass $M |
| uniaxialMaterial Elastic 1 $K |
| element truss 1 1 2 1 1 |
| # |
| # load pattern |
| # |
| if {$inputT == "Uniform" } { |
| set vel0 [expr -1.0*$U*$omega] |
| # set vel0 0.0 |
| puts "vel0 $vel0" |
| set factor [expr -1.0*$U*$omega*$omega] |
| set AccelSeries "Sine 0.0 $tFinish $period -factor $factor"; |
| pattern UniformExcitation 1 1 -accel $AccelSeries -vel0 $vel0 ; # create Unifform excitation |
| fix 1 1 |
| recorder Node -file DUniform.out -time -node 2 -dof 1 disp; # displacements of free nodes |
| recorder plot DUniform.out "DISP" 10 10 625 450 -dT [expr 10*0.01] -columns 1 2 -columns 1 3 |
| } else { |
| set DispSeries "Sine 0.0 $tFinish $period -factor $factor"; |
| pattern MultipleSupport 1 { |
| groundMotion 1 Plain -disp $DispSeries |
| imposedMotion 1 1 1 |
| }; |
| recorder Node -file DMultiSupport.out -time -node 1 2 -dof 1 disp; # displacements of free nodes |
| recorder plot DMultiSupport.out "DISP" 10 10 625 450 -dT [expr 10*0.01] -columns 1 2 -columns 1 3 |
| } |
| # |
| # Create analysis |
| # |
| # integrator HHT 1.0 |
| # integrator Newmark 0.25 0.5 |
| integrator Transient Trapezoidal |
| system ProfileSPD |
| algorithm Linear |
| constraints Transformation |
| analysis Transient |
| # |
| # Perform Analysis |
| # |
| set dt [expr $periodN/100] |
| set nSteps [expr int($tFinish/$dt)] |
| analyze $nSteps $dt |
| print node 1 2 |
| } |
| # |
| # set the parameters & use the procedure |
| # |
| set K 100.0 |
| set periodN 2.0 |
| set period 4.0 |
| set U 10.0 |
| set tFinish 6.0 |
| set inputT MultiSupport |
| doModel $K $periodN $U $period $tFinish $inputT |
| set inputT Uniform |
| doModel $K $periodN $U $period $tFinish $inputT |
| # |
| # create Matlab script: |
| # |
| set matlabF [open matlab.m w] |
| puts $matlabF "load DUniform.out" |
| puts $matlabF "load DMultiSupport.out" |
| puts $matlabF "a=DMultiSupport;" |
| puts $matlabF "b=DUniform;" |
| puts $matlabF "w=2*pi/$period" |
| puts $matlabF "wn=2*pi/$periodN" |
| puts $matlabF "U=$U" |
| puts $matlabF "exact=(U/(1-(w/wn)*(w/wn)))*sin(w*a(:,1)) - ((U*w/wn)/(1-(w/wn)*(w/wn)))*sin(wn*a(:,1));"; |
| puts $matlabF "plot(a(:,1),a(:,3),a(:,1),b(:,2),'*',a(:,1),a(:,2)+b(:,2),a(:,1),exact(:,1));" |
| puts $matlabF "grid on;" |
| puts $matlabF "legend(\"MultiSupport\",\"Uniform\",\"Ground+Uniform\",\"Exact\")"; |
| close $matlabF; |
| exit |
| /trunk/DEVELOPER/material/c/elasticPPC.c |
|---|
| New file |
| 0,0 → 1,179 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| /* |
| ** $Revision: 1.8 $ |
| ** $Date: 2009/01/13 01:00:52 $ |
| ** $Source: /usr/local/cvs/OpenSees/PACKAGES/NewMaterial/c/elasticPPC.c,v $ |
| ** Written: fmk |
| ** |
| ** Description: This file contains the implementation of elasticPP material |
| */ |
| #include "elementAPI.h" |
| #include <math.h> |
| #include <stdlib.h> |
| #include <stdio.h> |
| #define DBL_EPSILON 1e-18 |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| OPS_Export void |
| elasticPPC (matObj *thisObj, |
| modelState *model, |
| double *strain, |
| double *tang, |
| double *stress, |
| int *isw, |
| int *result) |
| { |
| *result = 0; |
| if (*isw == ISW_INIT) { |
| double dData[2]; |
| int iData[1]; |
| /* get the input data - tag? E? eyp? */ |
| int numData = 1; |
| OPS_GetIntInput(&numData, iData); |
| numData = 2; |
| OPS_GetDoubleInput(&numData, dData); |
| /* Allocate the element state */ |
| thisObj->tag = iData[0]; |
| thisObj->nParam = 2; /* E, eyp */ |
| thisObj->nState = 2; /* strain, ep */ |
| OPS_AllocateMaterial(thisObj); |
| thisObj->theParam[0] = dData[0]; |
| thisObj->theParam[1] = dData[1]; |
| *result = OPS_UNIAXIAL_MATERIAL_TYPE; |
| } else if (*isw == ISW_COMMIT) { |
| double trialStrain = thisObj->tState[0]; |
| double f, fYieldSurface; // yield function |
| double E = thisObj->theParam[0]; |
| double eyp = thisObj->theParam[1]; |
| double ep = thisObj->cState[1]; |
| double fyp = E*eyp; |
| double fyn = -fyp; |
| // compute trial stress |
| double sigtrial = E * ( trialStrain - ep ); |
| // evaluate yield function |
| if ( sigtrial >= 0.0 ) |
| f = sigtrial - fyp; |
| else |
| f = -sigtrial + fyn; |
| fYieldSurface = - E * DBL_EPSILON; |
| if ( f > fYieldSurface ) { |
| // plastic |
| if ( sigtrial > 0.0 ) { |
| ep += f / E; |
| } else { |
| ep -= f / E; |
| } |
| } |
| thisObj->cState[0] = trialStrain; |
| thisObj->cState[1] = ep; |
| } else if (*isw == ISW_REVERT) { |
| int i; |
| for (i=0; i<3; i++) |
| thisObj->tState[i] = thisObj->cState[i]; |
| } else if (*isw == ISW_REVERT_TO_START) { |
| int i; |
| for (i=0; i<2; i++) { |
| thisObj->cState[i] = 0.0; |
| thisObj->tState[i] = 0.0; |
| } |
| } else if (*isw == ISW_FORM_TANG_AND_RESID) { |
| double trialStrain = *strain; |
| double f; // yield function |
| double trialStress, trialTangent, fYieldSurface; |
| double E = thisObj->theParam[0]; |
| double eyp = thisObj->theParam[1]; |
| double ep = thisObj->cState[1]; |
| double fyp = E*eyp; |
| double fyn = -fyp; |
| // compute trial stress |
| double sigtrial = E * ( trialStrain - ep ); |
| // evaluate yield function |
| if ( sigtrial >= 0.0 ) |
| f = sigtrial - fyp; |
| else |
| f = -sigtrial + fyn; |
| 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; |
| } |
| thisObj->tState[0] = trialStrain; |
| *stress = trialStress; |
| *tang = trialTangent; |
| } |
| } |
| OPS_Export void |
| localInit() |
| { |
| OPS_Error("elasticPPC uniaxial material - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n", 1); |
| } |
| /trunk/DEVELOPER/material/c/elasticPlaneStressC.c |
|---|
| New file |
| 0,0 → 1,131 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| /* |
| ** $Revision$ |
| ** $Date$ |
| ** $Source$ |
| ** |
| ** Written: fmk |
| ** |
| ** Description: This file contains the implementation of elasticPP material |
| */ |
| #include "elementAPI.h" |
| #include <math.h> |
| #include <stdlib.h> |
| #include <stdio.h> |
| #define DBL_EPSILON 1e-18 |
| #ifdef _USRDLL |
| #include <windows.h> |
| #define OPS_Export _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| static int initFlag = 0; |
| OPS_Export void |
| elasticPlaneStressC(matObj *thisObj, |
| modelState *model, |
| double *strain, |
| double *tang, |
| double *stress, |
| int *isw, |
| int *result) |
| { |
| *result = 0; |
| if (*isw == ISW_INIT) { |
| if (initFlag == 0) { |
| OPS_Error("elasticPlaneStressC nd material - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n", 1); |
| initFlag = 1; |
| } |
| double dData[2]; |
| int iData[1]; |
| /* get the input data - tag? E? v? */ |
| int numData = 1; |
| OPS_GetIntInput(&numData, iData); |
| numData = 2; |
| OPS_GetDoubleInput(&numData, dData); |
| /* Allocate the element state */ |
| thisObj->tag = iData[0]; |
| thisObj->nParam = 2; /* E, v */ |
| thisObj->nState = 3; /* strain */ |
| OPS_AllocateMaterial(thisObj); |
| thisObj->theParam[0] = dData[0]; |
| thisObj->theParam[1] = dData[1]; |
| *result = OPS_PLANESTRESS_TYPE; |
| } else if (*isw == ISW_COMMIT) { |
| int i; |
| for (i=0; i<3; i++) |
| thisObj->cState[i] = thisObj->tState[i]; |
| } else if (*isw == ISW_REVERT) { |
| int i; |
| for (i=0; i<3; i++) |
| thisObj->tState[i] = thisObj->cState[i]; |
| } else if (*isw == ISW_REVERT_TO_START) { |
| int i; |
| for (i=0; i<3; i++) { |
| thisObj->cState[i] = 0.0; |
| thisObj->tState[i] = 0.0; |
| } |
| } else if (*isw == ISW_FORM_TANG_AND_RESID) { |
| int i; |
| for (i=0; i<3; i++) |
| thisObj->tState[i] = strain[i]; |
| // form tangent |
| double E = thisObj->theParam[0]; |
| double v = thisObj->theParam[1]; |
| double fact = E/(1-v*v); |
| tang[0] = fact; // tang[1,1] |
| tang[1] = fact*v; // tang[2,1] |
| tang[2] = 0.0; // tang[3,1] |
| tang[3] = fact*v; // tang[1,2] |
| tang[4] = fact; // tang[2,2] |
| tang[5] = 0.0; // tang[3,2] |
| tang[6] = 0.0; // tang[1,3] |
| tang[7] = 0.0; // tang[2,3] |
| tang[8] = fact*(1-v)/2.0; // tang[3,3] |
| stress[0] = tang[0]*strain[0] + tang[3]*strain[1]; |
| stress[1] = tang[1]*strain[0] + tang[4]*strain[1]; |
| stress[2] = 2.0 * tang[8] * strain[2]; |
| } |
| } |
| /trunk/DEVELOPER/material/c/Makefile |
|---|
| New file |
| 0,0 → 1,44 |
| include ../../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| all: $(LIBRARY_CPP) $(LIBRARY_C) $(LIBRARY_F) |
| LIBRARY_C = elasticPPc.so |
| OBJS_C = elasticPPc.o elasticPlaneStress.o |
| $(LIBRARY_C): $(OBJS_C) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY_CPP) $(OBJS_CPP) -lc |
| elasticPPC.o: elasticPPC.c |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall elasticPPC.c |
| elasticPlaneStressC.o: elasticPlaneStressC.c |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall elasticPlaneStress.c |
| else |
| LIBRARY_C = elasticPPC.dylib |
| all: $(OBJS) |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace elasticPPC.c $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o elasticPPC.dylib |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace elasticPlaneStressC.c $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o elasticPlaneStressC.dylib |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/material/c/example1.tcl |
|---|
| New file |
| 0,0 → 1,60 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.1 $ |
| # $Date: 2008/12/01 23:30:00 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewMaterial/c/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial elasticPPC 1 3000 0.03 |
| uniaxialMaterial elasticPPC 2 3000 0.001 |
| #uniaxialMaterial Elastic 1 3000.0 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element truss 1 1 4 10.0 1 |
| element truss 2 2 4 5.0 2 |
| element truss 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system BandSPD |
| constraints Plain |
| integrator LoadControl 1.0 1 1.0 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| exit |
| /trunk/DEVELOPER/material/c/example2.tcl |
|---|
| New file |
| 0,0 → 1,94 |
| # OpenSees Example: Simple supported beam modeled with 2D solid elements |
| # Units: kips, in, sec |
| # ---------------------------- |
| wipe; # clear opensees model |
| # create data directory |
| file mkdir Data |
| #----------------------------- |
| # Define the model |
| # ---------------------------- |
| # Create ModelBuilder with 2 dimensions and 2 DOF/node |
| model BasicBuilder -ndm 2 -ndf 2 |
| # create the material |
| nDMaterial elasticPlaneStressC 1 30000 0.25 |
| # set type of quadrilateral element (uncomment one of the three #options) |
| set Quad quad |
| #set Quad bbarQuad |
| #set Quad enhancedQuad |
| # set up the arguments for the three considered elements |
| if {$Quad == "enhancedQuad" } { |
| set eleArgs "PlaneStrain2D 1" |
| } |
| if {$Quad == "quad" } { |
| set eleArgs "1 PlaneStrain2D 1" |
| } |
| if {$Quad == "bbarQuad" } { |
| set eleArgs "1" |
| } |
| #define nodes |
| node 1 0.0 0.0 |
| node 2 100.0 0.0 |
| node 3 100.0 100.0 |
| node 4 0.0 100.0 |
| #define element |
| element quad 1 1 2 3 4 1.0 "PlaneStress" 1 |
| # define boundary conditions |
| fix 1 1 1 |
| fix 4 1 1 |
| # define the recorder |
| #--------------------- |
| recorder Node -file Data/RCUniaxialCompression.out -time -node 2 -dof disp |
| # define load pattern |
| #--------------------- |
| pattern Plain 1 Linear { |
| load 2 1.0 0.0 |
| load 3 1.0 0.0 |
| } |
| # -------------------------------------------------------------------- |
| # Start of static analysis (creation of the analysis & analysis itself) |
| # -------------------------------------------------------------------- |
| # Load control with variable load steps |
| # init Jd min max |
| integrator LoadControl 1.0 1 1.0 100.0 |
| # Convergence test |
| # tolerance maxIter displayCode |
| test EnergyIncr 1.0e-12 10 0 |
| # Solution algorithm |
| algorithm Newton |
| # DOF numberer |
| numberer RCM |
| # Cosntraint handler |
| constraints Plain |
| # System of equations solver |
| system ProfileSPD |
| # Type of analysis analysis |
| analysis Static |
| # Perform the analysis |
| analyze 10 |
| # -------------------------- |
| # End of static analysis |
| # -------------------------- |
| print ele |
| /trunk/DEVELOPER/material/cpp/ElasticPPCPP.cpp |
|---|
| New file |
| 0,0 → 1,283 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision: 1.7 $ |
| // $Date: 2009/03/23 23:17:04 $ |
| // $Source: /usr/local/cvs/OpenSees/PACKAGES/NewMaterial/cpp/ElasticPPcpp.cpp,v $ |
| // Written: fmk |
| // |
| // Description: This file contains the class implementation for |
| // ElasticMaterial. |
| // |
| // What: "@(#) ElasticPPcpp.C, revA" |
| #include <elementAPI.h> |
| #include "ElasticPPcpp.h" |
| #include <Vector.h> |
| #include <Channel.h> |
| #include <math.h> |
| #include <float.h> |
| #ifdef _USRDLL |
| #define OPS_Export extern "C" _declspec(dllexport) |
| #elif _MACOSX |
| #define OPS_Export extern "C" __attribute__((visibility("default"))) |
| #else |
| #define OPS_Export extern "C" |
| #endif |
| OPS_Export void |
| localInit() |
| { |
| OPS_Error("ElasticPPcpp unaxial material - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n", 1); |
| } |
| OPS_Export void * |
| OPS_ElasticPPcpp() |
| { |
| // Pointer to a uniaxial material that will be returned |
| UniaxialMaterial *theMaterial = 0; |
| int iData[1]; |
| double dData[2]; |
| int numData; |
| 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; |
| } |
| theMaterial = new ElasticPPcpp(iData[0], dData[0], dData[1]); |
| if (theMaterial == 0) { |
| opserr << "WARNING could not create uniaxialMaterial of type ElasticPPCpp\n"; |
| return 0; |
| } |
| return theMaterial; |
| } |
| ElasticPPcpp::ElasticPPcpp(int tag, double e, double eyp) |
| :UniaxialMaterial(tag,MAT_TAG_ElasticPPcpp), |
| ezero(0.0), E(e), trialStrain(0.0), ep(0.0), |
| trialStress(0.0), trialTangent(E) |
| { |
| fyp = E*eyp; |
| fyn = -fyp; |
| } |
| ElasticPPcpp::ElasticPPcpp() |
| :UniaxialMaterial(0,MAT_TAG_ElasticPPcpp), |
| fyp(0.0), fyn(0.0), ezero(0.0), E(0.0), trialStrain(0.0), ep(0.0), |
| trialStress(0.0), trialTangent(0.0) |
| { |
| } |
| ElasticPPcpp::~ElasticPPcpp() |
| { |
| // does nothing |
| } |
| int |
| 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 ); |
| //sigtrial = E * trialStrain; |
| //sigtrial -= E * ezero; |
| //sigtrial -= E * 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::commitState(void) |
| { |
| 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 ) { |
| // plastic |
| if ( sigtrial > 0.0 ) { |
| ep += f / E; |
| } else { |
| ep -= f / E; |
| } |
| } |
| return 0; |
| } |
| int |
| ElasticPPcpp::revertToLastCommit(void) |
| { |
| return 0; |
| } |
| int |
| ElasticPPcpp::revertToStart(void) |
| { |
| ep = 0.0; |
| return 0; |
| } |
| UniaxialMaterial * |
| ElasticPPcpp::getCopy(void) |
| { |
| ElasticPPcpp *theCopy = |
| new ElasticPPcpp(this->getTag(),E,fyp/E); |
| theCopy->ep = this->ep; |
| return theCopy; |
| } |
| int |
| ElasticPPcpp::sendSelf(int cTag, Channel &theChannel) |
| { |
| int res = 0; |
| static Vector data(6); |
| data(0) = this->getTag(); |
| data(1) = ep; |
| data(2) = E; |
| data(3) = ezero; |
| data(4) = fyp; |
| data(5) = fyn; |
| res = theChannel.sendVector(this->getDbTag(), cTag, data); |
| if (res < 0) |
| opserr << "ElasticPPcpp::sendSelf() - failed to send data\n"; |
| return res; |
| } |
| int |
| ElasticPPcpp::recvSelf(int cTag, Channel &theChannel, |
| FEM_ObjectBroker &theBroker) |
| { |
| int res = 0; |
| static Vector data(6); |
| res = theChannel.recvVector(this->getDbTag(), cTag, data); |
| if (res < 0) |
| opserr << "ElasticPPcpp::recvSelf() - failed to recv data\n"; |
| else { |
| this->setTag(data(0)); |
| ep = data(1); |
| E = data(2); |
| ezero = data(3); |
| fyp = data(4); |
| fyn = data(5); |
| } |
| return res; |
| } |
| void |
| ElasticPPcpp::Print(OPS_Stream &s, int flag) |
| { |
| s << "ElasticPPcpp tag: " << this->getTag() << endln; |
| s << " E: " << E << endln; |
| s << " ep: " << ep << endln; |
| s << " stress: " << trialStress << " tangent: " << trialTangent << endln; |
| } |
| /trunk/DEVELOPER/material/cpp/ElasticPPCPP.h |
|---|
| New file |
| 0,0 → 1,84 |
| /* ****************************************************************** ** |
| ** OpenSees - Open System for Earthquake Engineering Simulation ** |
| ** Pacific Earthquake Engineering Research Center ** |
| ** ** |
| ** ** |
| ** (C) Copyright 1999, The Regents of the University of California ** |
| ** All Rights Reserved. ** |
| ** ** |
| ** Commercial use of this program without express permission of the ** |
| ** University of California, Berkeley, is strictly prohibited. See ** |
| ** file 'COPYRIGHT' in main directory for information on usage and ** |
| ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
| ** ** |
| ** Developed by: ** |
| ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
| ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
| ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
| ** ** |
| ** ****************************************************************** */ |
| // $Revision: 1.1 $ |
| // $Date: 2008/12/09 20:00:16 $ |
| // $Source: /usr/local/cvs/OpenSees/PACKAGES/NewMaterial/cpp/ElasticPPcpp.h,v $ |
| #ifndef ElasticPPcpp_h |
| #define ElasticPPcpp_h |
| // Written: fmk |
| // |
| // Description: This file contains the class definition for |
| // ElasticPPcpp. ElasticPPcpp provides the abstraction |
| // of an elastic perfectly plastic uniaxial material, |
| // |
| // What: "@(#) ElasticPPcpp.h, revA" |
| #include <UniaxialMaterial.h> |
| #define MAT_TAG_ElasticPPcpp 5002 |
| class ElasticPPcpp : public UniaxialMaterial |
| { |
| public: |
| ElasticPPcpp(int tag, double E, double eyp); |
| ElasticPPcpp(); |
| ~ElasticPPcpp(); |
| int setTrialStrain(double strain, double strainRate = 0.0); |
| double getStrain(void); |
| double getStress(void); |
| double getTangent(void); |
| double getInitialTangent(void) {return E;}; |
| 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); |
| protected: |
| private: |
| double fyp, fyn; // positive and negative yield stress |
| double ezero; // initial strain |
| double E; // elastic modulus |
| double trialStrain; // trial strain |
| double ep; // plastic strain at last commit |
| double trialStress; // current trial stress |
| double trialTangent; // current trial tangent |
| }; |
| #endif |
| /trunk/DEVELOPER/material/cpp/Makefile |
|---|
| New file |
| 0,0 → 1,43 |
| include ../../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| all: $(LIBRARY_CPP) $(LIBRARY_C) $(LIBRARY_F) |
| LIBRARY_CPP = ElasticPPcpp.so |
| OBJS_CPP = ElasticPPcpp.o |
| $(LIBRARY_CPP): $(OBJS_CPP) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY_CPP) $(OBJS_CPP) -lc |
| ElasticPPcpp.o: ElasticPPcpp.cpp |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall ElasticPPcpp.cpp |
| elasticPPC.o: elasticPPC.c |
| $(CC++) -fPIC $(INCLUDES) -g -c -Wall elasticPPC.c |
| else |
| LIBRARY_CPP = ElasticPPCPP.dylib |
| all: $(OBJS) |
| $(CC++) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace ElasticPPcpp.cpp $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY_CPP) |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/material/cpp/example1.tcl |
|---|
| New file |
| 0,0 → 1,58 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.2 $ |
| # $Date: 2008/12/09 19:59:46 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewMaterial/cpp/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial ElasticPPcpp 1 3000 0.001 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element truss 1 1 4 10.0 1 |
| element truss 2 2 4 5.0 1 |
| element truss 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system BandSPD |
| constraints Plain |
| integrator LoadControl 1.0 1 1.0 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| exit |
| /trunk/DEVELOPER/material/fortran/materialTypes.f |
|---|
| New file |
| 0,0 → 1,45 |
| module materialTypes |
| use iso_c_binding |
| implicit none |
| private |
| public modelState |
| type, bind(c) :: modelState |
| real (C_DOUBLE) time |
| real (C_DOUBLE) dt |
| end type modelState |
| c type modelState |
| c real *8::time; |
| c real *8::dt; |
| c end type modelState |
| public matObject |
| type, bind(C) :: matObject |
| integer(C_INT) ::tag; |
| integer(C_INT) ::nParam; |
| integer(C_INT) ::nState |
| type (c_ptr) :: theParam; |
| type (c_ptr) :: cState; |
| type (c_ptr) :: tState; |
| type(c_funptr) ::functPtr |
| type (c_ptr) :: matObjPtr; |
| end type matObject; |
| c real (C_DOUBLE) :: theParam; |
| c real (C_DOUBLE) :: cState; |
| c real (C_DOUBLE) :: tState; |
| c public matObject |
| c type :: matObject |
| c integer::tag |
| c integer::nParam |
| c integer::nState |
| c real *8, pointer ::theParam(:) |
| c real *8, pointer ::cState(:) |
| c real *8, pointer ::tState(:) |
| c type(c_funptr)::functPtr |
| c end type matObject |
| end module |
| /trunk/DEVELOPER/material/fortran/materialAPI.f |
|---|
| New file |
| 0,0 → 1,65 |
| module materialAPI |
| use iso_c_binding |
| use materialTypes |
| implicit none |
| integer ::ISW_COMMIT, ISW_INIT, ISW_REVERT_TO_START |
| integer ::ISW_REVERT, ISW_FORM_TANG_AND_RESID, ISW_FORM_MASS |
| real *8 ::DBL_EPSILON |
| PARAMETER (ISW_INIT = 0) |
| PARAMETER (ISW_COMMIT = 1) |
| PARAMETER (ISW_REVERT = 2) |
| PARAMETER (ISW_FORM_TANG_AND_RESID = 3) |
| PARAMETER (ISW_FORM_MASS = 4) |
| PARAMETER (ISW_REVERT_TO_START = 5) |
| PARAMETER (DBL_EPSILON = 1.0e-21) |
| c PUBLIC OPS_GetIntInput |
| c interface |
| c function OPS_GetIntInput(numData, iData) |
| c integer :: OPS_GetIntInput |
| c integer :: numData |
| c integer :: iData(:) |
| c end function OPS_GetIntInput |
| c end interface |
| PUBLIC OPS_GetIntInput |
| interface |
| function OPS_GetIntInput(numData, iData) |
| integer :: OPS_GetIntInput |
| integer :: numData |
| integer :: iData(*) |
| end function OPS_GetIntInput |
| end interface |
| PUBLIC OPS_GetDoubleInput |
| interface |
| function OPS_GetDoubleInput(numData, dData) |
| integer :: OPS_GetDoubleInput |
| integer :: numData |
| real *8 :: dData(*) |
| end function OPS_GetDoubleInput |
| end interface |
| PUBLIC OPS_AllocateMaterial |
| interface |
| function OPS_AllocateMaterial(mat) |
| use materialtypes |
| integer :: OPS_AllocateMaterial |
| type(matObject) :: mat |
| end function OPS_AllocateMaterial |
| end interface |
| PUBLIC OPS_Error |
| interface |
| function OPS_Error(msg) |
| integer :: OPS_Error |
| character, dimension(*) :: msg |
| end function OPS_Error |
| end interface |
| c contains |
| end module materialAPI |
| /trunk/DEVELOPER/material/fortran/Makefile |
|---|
| New file |
| 0,0 → 1,47 |
| include ../../Makefile.def |
| ifeq ($(OS_FLAG), -D_LINUX) |
| all: $(LIBRARY_F) |
| LIBRARY_F = elasticPPf.so |
| OBJS_F = elasticPPF.o |
| $(LIBRARY_F): $(OBJS_F) |
| $(CC++) -shared -Wl,-soname,$(LIBRARY) -o $(LIBRARY) elasticPPF.o -lc |
| else |
| LIBRARY_F = elasticPPf.dylib |
| OTHER_OBJS = dynamicLoadFunctions.o |
| OBJS = materialtypes.mod materialapi.mod |
| all: $(OBJS) |
| $(FC) $(OS_FLAG) -dynamiclib $(INCLUDES) -Wl,-undefined,suppress,-flat_namespace elasticPPf.f $(OUTSIDE_OBJS) -current_version 1.0 -compatibility_version 1.0 -fvisibility=hidden -o $(LIBRARY_F) |
| materialapi.mod: materialAPI.f materialtypes.mod |
| $(FC) -c materialAPI.f |
| materialtypes.mod: materialTypes.f |
| $(FC) -c materialTypes.f |
| endif |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core |
| wipe: spotless |
| @$(RM) $(RMFLAGS) $(PROGRAM) fake core $(LIBRARY) |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |
| /trunk/DEVELOPER/material/fortran/example1.tcl |
|---|
| New file |
| 0,0 → 1,59 |
| # written: fmk |
| # |
| # purpose: example1 in the manual with existing truss element |
| # |
| # $Revision: 1.2 $ |
| # $Date: 2009/01/12 20:37:57 $ |
| # $Source: /usr/local/cvs/OpenSees/PACKAGES/NewMaterial/fortran/example1.tcl,v $ |
| #create the ModelBuilder object |
| model basic -ndm 2 -ndf 2 |
| # build the model |
| # add nodes - command: node nodeId xCrd yCrd |
| node 1 0.0 0.0 |
| node 2 144.0 0.0 |
| node 3 168.0 0.0 |
| node 4 72.0 96.0 |
| # add material - command: material <matType> matID <matArgs> |
| uniaxialMaterial elasticppf 1 3000 0.03 |
| uniaxialMaterial elasticppf 2 3000 0.001 |
| #uniaxialMaterial Elastic 1 3000.0 |
| # add truss elements - command: truss trussID node1 node2 A matID |
| element truss 1 1 4 10.0 1 |
| element truss 2 2 4 5.0 2 |
| element truss 3 3 4 5.0 1 |
| # set the boundary conditions - command: fix nodeID xResrnt? yRestrnt? |
| fix 1 1 1 |
| fix 2 1 1 |
| fix 3 1 1 |
| pattern Plain 1 Linear { |
| # apply the load - command: load nodeID xForce yForce |
| load 4 100 -50 |
| } |
| # build the components for the analysis object |
| system BandSPD |
| constraints Plain |
| integrator LoadControl 1.0 1 1.0 1.0 |
| algorithm Linear |
| numberer RCM |
| # create the analysis object |
| analysis Static |
| # create a Recorder object for the nodal displacements at node 4 |
| recorder Node -file Example.out -load -nodes 4 -dof 1 2 disp |
| # perform the analysis |
| analyze 1 |
| # print the results at node and at all elements |
| print node 4 |
| print ele |
| /trunk/DEVELOPER/material/fortran/elasticPPf.f |
|---|
| New file |
| 0,0 → 1,168 |
| SUBROUTINE ELASTICPPF(matObj,model,strain,tang,stress,isw,error) |
| !DEC$ IF DEFINED (_DLL) |
| !DEC$ ATTRIBUTES DLLEXPORT :: ELASTICPPF |
| !DEC$ END IF |
| use materialTypes |
| use materialAPI |
| implicit none |
| type(matObject)::matObj |
| type(modelState)::model |
| real *8::strain |
| real *8::tang |
| real *8::stress |
| integer::isw |
| integer::error; |
| real *8, pointer::theParam(:) |
| real *8, pointer::cState(:) |
| real *8, pointer::tState(:) |
| real *8 E, eyp, ep, fyp, fyn, f, fYieldSurface |
| real *8 sigtrial, trialStrain, trialStress, trialTangent |
| integer, target :: iData(2) |
| integer, pointer ::iPtr(:) |
| integer numData, err |
| real *8, target :: dData(2) |
| real *8, pointer:: dPtr(:) |
| c outside functions called |
| c integer OPS_GetIntInput, OPS_GetDoubleInput, OPS_AllocateMaterial |
| IF (isw.eq.ISW_INIT) THEN |
| c get the input data - tag? E? eyp? |
| numData = 1 |
| iPtr=>iData; |
| err = OPS_GetIntInput(numData, iPtr) |
| numData = 2 |
| dPtr=>dData; |
| err = OPS_GetDoubleInput(numData, dPtr) |
| c Allocate the element state |
| matObj%tag = idata(1) |
| matObj%nparam = 2 |
| matObj%nstate = 2 |
| err = OPS_AllocateMaterial(matObj) |
| call c_f_pointer(matObj%theParam, theParam, [2]); |
| c Initialize the element properties |
| theParam(1) = dData(1); |
| theParam(2) = dData(2); |
| ELSE |
| call c_f_pointer(matObj%theParam, theParam, [2]); |
| call c_f_pointer(matObj%cState, cState, [2]); |
| call c_f_pointer(matObj%tState, tState, [2]); |
| IF (isw == ISW_COMMIT) THEN |
| trialStrain = tState(1) |
| E = theParam(1) |
| eyp = theParam(2) |
| ep = cState(2) |
| fyp = E*eyp |
| fyn = -fyp |
| c compute trial stress |
| sigtrial = E * ( trialStrain - ep ); |
| c evaluate yield function |
| IF ( sigtrial >= 0.0 ) THEN |
| f = sigtrial - fyp |
| ELSE |
| f = -sigtrial + fyn |
| END IF |
| fYieldSurface = - E * DBL_EPSILON; |
| IF ( f > fYieldSurface ) THEN |
| c plastic |
| IF ( sigtrial > 0.0 ) THEN |
| ep = ep + f / E |
| ELSE |
| ep = ep - f / E |
| END IF |
| END IF |
| cState(1) = trialStrain; |
| cState(2) = ep; |
| ELSE IF (isw == ISW_REVERT_TO_START) THEN |
| cState(1) = 0.0; |
| cState(2) = 0.0; |
| tState(1) = 0.0; |
| tState(2) = 0.0; |
| ELSE IF (isw == ISW_FORM_TANG_AND_RESID) THEN |
| trialStrain = strain; |
| E = theParam(1); |
| eyp = theParam(2); |
| ep = cState(2); |
| fyp = E*eyp; |
| fyn = -fyp; |
| c compute trial stress |
| sigtrial = E * ( trialStrain - ep ); |
| c evaluate yield function |
| IF ( sigtrial >= 0.0 ) THEN |
| f = sigtrial - fyp; |
| ELSE |
| f = -sigtrial + fyn; |
| END IF |
| fYieldSurface = - E * DBL_EPSILON; |
| IF ( f <= fYieldSurface ) THEN |
| trialStress = sigtrial; |
| trialTangent = E; |
| ELSE |
| c plastic |
| IF ( sigtrial > 0.0 ) THEN |
| trialStress = fyp; |
| ELSE |
| trialStress = fyn; |
| END IF |
| trialTangent = 0.0; |
| END IF |
| tState(1) = trialStrain; |
| stress = trialStress; |
| tang = trialTangent; |
| END IF |
| END IF |
| c return error code |
| error = 0 |
| END SUBROUTINE elasticPPf |
| SUBROUTINE LOCALINIT() |
| !DEC$ IF DEFINED (_DLL) |
| !DEC$ ATTRIBUTES DLLEXPORT :: LOCALINIT |
| !DEC$ END IF |
| use materialAPI |
| implicit none |
| integer::error; |
| character *60:: msg; |
| msg = 'elasticPPf - Written by fmk UC Berkeley Copyright 2008' |
| error = OPS_Error(msg); |
| END SUBROUTINE localInit |
| /trunk/DEVELOPER/material/Makefile |
|---|
| New file |
| 0,0 → 1,24 |
| include ../Makefile.def |
| all: |
| echo @pwd; |
| @$(CD) cpp; $(MAKE); |
| @$(CD) c; $(MAKE); |
| @$(CD) fortran; $(MAKE); |
| # Miscellaneous |
| tidy: |
| @$(RM) $(RMFLAGS) Makefile.bak *~ #*# core |
| clean: tidy |
| @$(CD) ./cpp; $(MAKE) wipe; |
| @$(CD) ./c; $(MAKE) wipe; |
| @$(CD) ./fortran; $(MAKE) wipe; |
| @$(RM) $(RMFLAGS) $(OBJS) *.o core *.out *.so *.dylib |
| spotless: clean |
| wipe: spotless |
| # DO NOT DELETE THIS LINE -- make depend depends on it. |