diff --git a/DEVELOPER/core/classTags.h b/DEVELOPER/core/classTags.h index d356c6345f..65de56e9ce 100644 --- a/DEVELOPER/core/classTags.h +++ b/DEVELOPER/core/classTags.h @@ -810,6 +810,7 @@ #define HANDLER_TAG_PenaltyConstraintHandler 3 #define HANDLER_TAG_TransformationConstraintHandler 4 #define HANDLER_TAG_PenaltyHandlerNoHomoSPMultipliers 5 +#define HANDLER_TAG_AutoConstraintHandler 6 #define NUMBERER_TAG_DOF_Numberer 1 #define NUMBERER_TAG_PlainNumberer 2 diff --git a/SRC/Makefile b/SRC/Makefile index 37583d6864..9d6ba4aa3c 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -1489,6 +1489,7 @@ ANALYSIS_LIBS = $(FE)/analysis/analysis/Analysis.o \ $(FE)/analysis/handler/PenaltyConstraintHandler.o \ $(FE)/analysis/handler/LagrangeConstraintHandler.o \ $(FE)/analysis/handler/TransformationConstraintHandler.o \ + $(FE)/analysis/handler/AutoConstraintHandler.o \ $(FE)/analysis/numberer/DOF_Numberer.o \ $(FE)/analysis/numberer/PlainNumberer.o \ $(FE)/analysis/numberer/ParallelNumberer.o \ diff --git a/SRC/actor/objectBroker/FEM_ObjectBrokerAllClasses.cpp b/SRC/actor/objectBroker/FEM_ObjectBrokerAllClasses.cpp index e227a81fdd..5a55ec31ce 100644 --- a/SRC/actor/objectBroker/FEM_ObjectBrokerAllClasses.cpp +++ b/SRC/actor/objectBroker/FEM_ObjectBrokerAllClasses.cpp @@ -587,6 +587,7 @@ #include "PenaltyConstraintHandler.h" #include "LagrangeConstraintHandler.h" #include "TransformationConstraintHandler.h" +#include "AutoConstraintHandler.h" // dof numberer header files #include "DOF_Numberer.h" @@ -2704,7 +2705,10 @@ FEM_ObjectBrokerAllClasses::getNewConstraintHandler(int classTag) case HANDLER_TAG_TransformationConstraintHandler: return new TransformationConstraintHandler(); - + + case HANDLER_TAG_AutoConstraintHandler: + return new AutoConstraintHandler(); + default: opserr << "FEM_ObjectBrokerAllClasses::getNewConstraintHandler - "; opserr << " - no ConstraintHandler type exists for class tag "; diff --git a/SRC/analysis/dof_grp/LagrangeDOF_Group.cpp b/SRC/analysis/dof_grp/LagrangeDOF_Group.cpp index 60b99ea464..a85d03cbff 100644 --- a/SRC/analysis/dof_grp/LagrangeDOF_Group.cpp +++ b/SRC/analysis/dof_grp/LagrangeDOF_Group.cpp @@ -42,14 +42,23 @@ #include #include +//#define LAG_DOF_VERBOSE + LagrangeDOF_Group::LagrangeDOF_Group(int tag, SP_Constraint &spPtr) :DOF_Group(tag, 1) +, m_lagrange_variable(1) { } LagrangeDOF_Group::LagrangeDOF_Group(int tag, MP_Constraint &mpPtr) -:DOF_Group(tag, (mpPtr.getRetainedDOFs()).Size()) +// changed by M.Petracca: 2024. should be the constrained dof size. +// now it's the same because all MP constraints in opensees are one-to-one. +// but in the future we may implement generic constraints of the form: +// CDOF = a1*RDOF1 + a2*RDOF2 + ... + an*RDOFn + beta +//:DOF_Group(tag, (mpPtr.getRetainedDOFs()).Size()) +:DOF_Group(tag, (mpPtr.getConstrainedDOFs()).Size()) +, m_lagrange_variable(mpPtr.getConstrainedDOFs().Size()) { } @@ -98,7 +107,17 @@ LagrangeDOF_Group::getUnbalance(Integrator *theIntegrator) void LagrangeDOF_Group::setNodeDisp(const Vector &u) { - return; + m_lagrange_variable.Zero(); + const auto& ids = this->getID(); + for (int i = 0; i < this->getNumDOF(); i++) { + int loc = ids(i); + if (loc >= 0) + m_lagrange_variable(i) = u(loc); + } +#ifdef LAG_DOF_VERBOSE + opserr << "LAG DOF: set U (total) = " << u; + opserr << " -> LM = " << m_lagrange_variable; +#endif // LAG_DOF_VERBOSE } @@ -132,7 +151,16 @@ LagrangeDOF_Group::setNodeAccel(const Vector &udotdot) void LagrangeDOF_Group::incrNodeDisp(const Vector &u) { - return; + const auto& ids = this->getID(); + for (int i = 0; i < this->getNumDOF(); i++) { + int loc = ids(i); + if (loc >= 0) + m_lagrange_variable(i) += u(loc); + } +#ifdef LAG_DOF_VERBOSE + opserr << "LAG DOF: set U (increment) = " << u; + opserr << " -> LM = " << m_lagrange_variable; +#endif // LAG_DOF_VERBOSE } @@ -162,8 +190,10 @@ LagrangeDOF_Group::incrNodeAccel(const Vector &udotdot) const Vector & LagrangeDOF_Group::getCommittedDisp(void) { - unbalance->Zero(); - return *unbalance; + // note: this is actually the trial one. but this method is only + // called by triansient integrators during the domainChanged method + // (trial and committed should be the same) + return m_lagrange_variable; } const Vector & @@ -178,6 +208,23 @@ LagrangeDOF_Group::getCommittedAccel(void) { unbalance->Zero(); return *unbalance; +} + +const Vector& LagrangeDOF_Group::getTrialDisp() +{ + return m_lagrange_variable; +} + +const Vector& LagrangeDOF_Group::getTrialVel() +{ + unbalance->Zero(); + return *unbalance; +} + +const Vector& LagrangeDOF_Group::getTrialAccel() +{ + unbalance->Zero(); + return *unbalance; } void diff --git a/SRC/analysis/dof_grp/LagrangeDOF_Group.h b/SRC/analysis/dof_grp/LagrangeDOF_Group.h index 44de52f7f7..57b90115c9 100644 --- a/SRC/analysis/dof_grp/LagrangeDOF_Group.h +++ b/SRC/analysis/dof_grp/LagrangeDOF_Group.h @@ -39,6 +39,7 @@ // What: "@(#) LagrangeDOF_Group.h, revA" #include +#include class SP_Constraint; class MP_Constraint; @@ -59,6 +60,9 @@ class LagrangeDOF_Group: public DOF_Group virtual const Vector &getCommittedDisp(void); virtual const Vector &getCommittedVel(void); virtual const Vector &getCommittedAccel(void); + virtual const Vector &getTrialDisp(); + virtual const Vector &getTrialVel(); + virtual const Vector &getTrialAccel(); // methods to update the trial response at the nodes virtual void setNodeDisp(const Vector &u); @@ -83,6 +87,9 @@ class LagrangeDOF_Group: public DOF_Group protected: private: + // we don't have a physical Node, so we need a persistent storage + // for the lagrange multipliers so that the lagrange FE can correctly compute the residual + Vector m_lagrange_variable; }; diff --git a/SRC/analysis/dof_grp/TransformationDOF_Group.cpp b/SRC/analysis/dof_grp/TransformationDOF_Group.cpp index ab0ea83251..4ca81bc684 100644 --- a/SRC/analysis/dof_grp/TransformationDOF_Group.cpp +++ b/SRC/analysis/dof_grp/TransformationDOF_Group.cpp @@ -169,6 +169,11 @@ TransformationDOF_Group::TransformationDOF_Group(int tag, Node *node, exit(-1); } } + +#ifdef TRANSF_INCREMENTAL_MP + modTotalDisp.resize(modNumDOF); + modTotalDisp = getTrialDisp(); +#endif // TRANSF_INCREMENTAL_MP numTransDOFs++; theHandler = theTHandler; @@ -227,6 +232,11 @@ TransformationDOF_Group::TransformationDOF_Group(int tag, } } +#ifdef TRANSF_INCREMENTAL_MP + modTotalDisp.resize(modNumDOF); + modTotalDisp = getTrialDisp(); +#endif // TRANSF_INCREMENTAL_MP + numTransDOFs++; theHandler = theTHandler; } @@ -469,6 +479,12 @@ TransformationDOF_Group::getCommittedAccel(void) void TransformationDOF_Group::setNodeDisp(const Vector &u) { +#ifdef TRANSF_INCREMENTAL_MP + // save the previous mod trial here + static Vector modTrialDispOld; + modTrialDispOld = modTotalDisp; // at previous iteration +#endif // TRANSF_INCREMENTAL_MP + // call base class method and return if no MP_Constraint if (theMP == 0) { this->DOF_Group::setNodeDisp(u); @@ -499,6 +515,31 @@ TransformationDOF_Group::setNodeDisp(const Vector &u) } } +#ifdef TRANSF_INCREMENTAL_MP + modTotalDisp = *modUnbalance; // save it for next iteration +#endif // TRANSF_INCREMENTAL_MP + + // at this point the modUnbalance contains the reduced total displacement. + // remove the trial one to obtain the increment, so that we transform only the increment +#ifdef TRANSF_INCREMENTAL_MP + modUnbalance->addVector(1.0, modTrialDispOld, -1.0); +#ifdef TRANSF_INCREMENTAL_MP_DEBUG + opserr << " N = " << myNode->getTag() << "\n"; + opserr << " solut: " << u; + opserr << " oldtr: " << modTrialDispOld; + opserr << " trial: " << myNode->getTrialDisp(); + opserr << " commi: " << myNode->getDisp(); + opserr << " incre: " << myNode->getIncrDisp(); + opserr << " delta: " << myNode->getIncrDeltaDisp(); + Domain* dom = myNode->getDomain(); + Node* ret = dom->getNode(theMP->getNodeRetained()); + opserr << " R VEL: " << ret->getTrialVel(); + opserr << " C VEL: " << ret->getTrialVel(); + opserr << " R ACC: " << ret->getTrialAccel(); + opserr << " C ACC: " << ret->getTrialAccel(); +#endif // TRANSF_INCREMENTAL_MP_DEBUG +#endif // TRANSF_INCREMENTAL_MP + Matrix *T = this->getT(); // *unbalance = (*T) * (*modUnbalance); unbalance->addMatrixVector(0.0, *T, *modUnbalance, 1.0); @@ -508,9 +549,18 @@ TransformationDOF_Group::setNodeDisp(const Vector &u) int numDOF = myNode->getNumberDOF(); for (int i=0; iincrTrialDisp(*unbalance); +#else myNode->setTrialDisp(*unbalance); +#endif // #ifdef TRANSF_INCREMENTAL_MP } void @@ -628,6 +678,10 @@ TransformationDOF_Group::incrNodeDisp(const Vector &u) (*modUnbalance)(i) = 0.0; } +#ifdef TRANSF_INCREMENTAL_MP + modTotalDisp.addVector(1.0, *modUnbalance, 1.0); // accumulate it for next iteration +#endif // TRANSF_INCREMENTAL_MP + Matrix *T = this->getT(); // *unbalance = (*T) * (*modUnbalance); unbalance->addMatrixVector(0.0, *T, *modUnbalance, 1.0); @@ -985,7 +1039,16 @@ TransformationDOF_Group::enforceSPs(int doMP) for (int i=0; igetValue(); - myNode->setTrialDisp(value, i); +#ifdef TRANSF_INCREMENTAL_SP + // include the initial value for staged analyses. + // note: do it only here. No need to do it when doMP == 0 + // because it will be called after doMP == 1 and the value + // has already been set to the retained node + double initial_value = theSPs[i]->getInitialValue(); + myNode->setTrialDisp(value + initial_value, i); +#else + myNode->setTrialDisp(value, i); +#endif // TRANSF_INCREMENTAL_SP } } diff --git a/SRC/analysis/dof_grp/TransformationDOF_Group.h b/SRC/analysis/dof_grp/TransformationDOF_Group.h index a40ed5329c..891db366ab 100644 --- a/SRC/analysis/dof_grp/TransformationDOF_Group.h +++ b/SRC/analysis/dof_grp/TransformationDOF_Group.h @@ -40,6 +40,15 @@ #include +// M.Petracca 2024. Unified approach to constraints +#define TRANSF_INCREMENTAL_SP +#define TRANSF_INCREMENTAL_MP +//#define TRANSF_INCREMENTAL_MP_DEBUG + +#ifdef TRANSF_INCREMENTAL_MP +#include +#endif // TRANSF_INCREMENTAL_MP + class MP_Constraint; class SP_Constraint; class TransformationConstraintHandler; @@ -129,6 +138,14 @@ class TransformationDOF_Group: public DOF_Group static Vector **modVectors; // array of pointers to class widde vectors static int numTransDOFs; // number of objects static TransformationConstraintHandler *theHandler; + +#ifdef TRANSF_INCREMENTAL_MP + // used to store locally the total displacement as we cannot rely + // on getting it from the retained node when processing the constrained node: + // the retained node may have been processed before + Vector modTotalDisp; +#endif // TRANSF_INCREMENTAL_MP + }; #endif diff --git a/SRC/analysis/fe_ele/lagrange/LagrangeMP_FE.cpp b/SRC/analysis/fe_ele/lagrange/LagrangeMP_FE.cpp index af48ed1f1c..86acfa2ce3 100644 --- a/SRC/analysis/fe_ele/lagrange/LagrangeMP_FE.cpp +++ b/SRC/analysis/fe_ele/lagrange/LagrangeMP_FE.cpp @@ -46,11 +46,17 @@ #include #include +// Note 1: +// changed by M.Petracca: 2024. should be the constrained dof size. +// now it's the same because all MP constraints in opensees are one-to-one. +// but in the future we may implement generic constraints of the form: +// CDOF = a1*RDOF1 + a2*RDOF2 + ... + an*RDOFn + beta + LagrangeMP_FE::LagrangeMP_FE(int tag, Domain &theDomain, MP_Constraint &TheMP, DOF_Group &theGroup, double Alpha) :FE_Element(tag, 3,(tag, TheMP.getConstrainedDOFs()).Size()+ (TheMP.getRetainedDOFs()).Size() + - (TheMP.getRetainedDOFs()).Size()), + (TheMP.getConstrainedDOFs()).Size()), // *see note 1 alpha(Alpha), theMP(&TheMP), theConstrainedNode(0), theRetainedNode(0), theDofGroup(&theGroup), tang(0), resid(0) @@ -104,6 +110,12 @@ LagrangeMP_FE::LagrangeMP_FE(int tag, Domain &theDomain, MP_Constraint &TheMP, exit(-1); } + if (theDofGroup->getID().Size() != theConstrainedNodesDOFs->getID().Size()) { + opserr << "WARNING LagrangeMP_FE::LagrangeMP_FE()"; + opserr << " - ConstrainedDOFs size != Lagrange size\n"; + exit(-1); + } + myDOF_Groups(0) = theConstrainedNodesDOFs->getTag(); myDOF_Groups(1) = theRetainedNodesDOFs->getTag(); myDOF_Groups(2) = theDofGroup->getTag(); @@ -224,6 +236,54 @@ LagrangeMP_FE::getTangent(Integrator *theNewIntegrator) const Vector & LagrangeMP_FE::getResidual(Integrator *theNewIntegrator) { + // get the solution vector [Uc Ur lambda] + static Vector UU; + const ID& id1 = theMP->getConstrainedDOFs(); + const ID& id2 = theMP->getRetainedDOFs(); + const ID& id3 = theDofGroup->getID(); + int size = id1.Size() + id2.Size() + id3.Size(); + UU.resize(size); + const Vector& Uc = theConstrainedNode->getTrialDisp(); + const Vector& Ur = theRetainedNode->getTrialDisp(); + const Vector& Uc0 = theMP->getConstrainedDOFsInitialDisplacement(); + const Vector& Ur0 = theMP->getRetainedDOFsInitialDisplacement(); + const Vector& lambda = theDofGroup->getTrialDisp(); + for (int i = 0; i < id1.Size(); ++i) { + int cdof = id1(i); + if (cdof < 0 || cdof >= Uc.Size()) { + opserr << "PenaltyMP_FE::getResidual FATAL Error: Constrained DOF " << cdof << " out of bounds [0-" << Uc.Size() << "]\n"; + exit(-1); + } + UU(i) = Uc(cdof) - Uc0(i); + } + for (int i = 0; i < id2.Size(); ++i) { + int rdof = id2(i); + if (rdof < 0 || rdof >= Ur.Size()) { + opserr << "PenaltyMP_FE::getResidual FATAL Error: Retained DOF " << rdof << " out of bounds [0-" << Ur.Size() << "]\n"; + exit(-1); + } + UU(i + id1.Size()) = Ur(rdof) - Ur0(i); + } + for (int i = 0; i < id3.Size(); ++i) { + UU(i + id1.Size() + id2.Size()) = lambda(i); + } + + /* + R = -C*U + G + .R = generalized residual vector + .C = constraint matrix + .U = generalized solution vector (displacement, lagrange multipliers) + .G = constrain values for non-homogeneous MP constraints (not available now) + | Ru | | 0 A | | u | | 0 | + | | = -| |*| | + | | + | Rl | | A 0 | | l | | g | + */ + + // compute residual + const Matrix& KK = getTangent(theNewIntegrator); + resid->addMatrixVector(0.0, KK, UU, -1.0); + + // done return *resid; } diff --git a/SRC/analysis/fe_ele/lagrange/LagrangeSP_FE.cpp b/SRC/analysis/fe_ele/lagrange/LagrangeSP_FE.cpp index 8538feb0dc..50e4bacab7 100644 --- a/SRC/analysis/fe_ele/lagrange/LagrangeSP_FE.cpp +++ b/SRC/analysis/fe_ele/lagrange/LagrangeSP_FE.cpp @@ -140,18 +140,36 @@ const Vector & LagrangeSP_FE::getResidual(Integrator *theNewIntegrator) { double constraint = theSP->getValue(); + double initialValue = theSP->getInitialValue(); int constrainedDOF = theSP->getDOF_Number(); const Vector &nodeDisp = theNode->getTrialDisp(); + const Vector& lambda = theDofGroup->getTrialDisp(); if (constrainedDOF < 0 || constrainedDOF >= nodeDisp.Size()) { - opserr << "LagrangeSP_FE::formResidual() -"; - opserr << " constrained DOF " << constrainedDOF << " outside range\n"; - (*resid)(1) = 0; + opserr << "LagrangeSP_FE::getResidual() -"; + opserr << " constrained DOF " << constrainedDOF << " outside range\n"; + resid->Zero(); + return *resid; + } + if (lambda.Size() != 1) { + opserr << "LagrangeSP_FE::getResidual() -"; + opserr << " Lambda.Size() = " << lambda.Size() << " != 1\n"; + resid->Zero(); + return *resid; } - (*resid)(1) = alpha *(constraint - nodeDisp(constrainedDOF)); -// opserr << "LagrangeSP_FE::getResidual() " << constraint << " " << nodeDisp(constrainedDOF) << " " << constrainedDOF << nodeDisp; -// opserr << "LagrangeSP_FE::getResidual() " << *resid << this->getID(); + /* + R = -C*U + G + .R = generalized residual vector + .C = constraint matrix + .U = generalized solution vector (displacement, lagrange multipliers) + .G = imposed displacement values + | Ru | | 0 A | | u | | 0 | + | | = -| |*| | + | | + | Rl | | A 0 | | l | | g | + */ + (*resid)(0) = alpha * (-lambda(0)); + (*resid)(1) = alpha *(constraint - (nodeDisp(constrainedDOF) - initialValue)); return *resid; } diff --git a/SRC/analysis/fe_ele/penalty/PenaltyMP_FE.cpp b/SRC/analysis/fe_ele/penalty/PenaltyMP_FE.cpp index 82724a16b4..db3ba8f719 100644 --- a/SRC/analysis/fe_ele/penalty/PenaltyMP_FE.cpp +++ b/SRC/analysis/fe_ele/penalty/PenaltyMP_FE.cpp @@ -210,6 +210,39 @@ const Vector & PenaltyMP_FE::getResidual(Integrator *theNewIntegrator) { // zero residual, CD = 0 + + // get the solution vector [Uc Ur] + static Vector UU; + const ID& id1 = theMP->getConstrainedDOFs(); + const ID& id2 = theMP->getRetainedDOFs(); + int size = id1.Size() + id2.Size(); + UU.resize(size); + const Vector& Uc = theConstrainedNode->getTrialDisp(); + const Vector& Ur = theRetainedNode->getTrialDisp(); + const Vector& Uc0 = theMP->getConstrainedDOFsInitialDisplacement(); + const Vector& Ur0 = theMP->getRetainedDOFsInitialDisplacement(); + for (int i = 0; i < id1.Size(); ++i) { + int cdof = id1(i); + if (cdof < 0 || cdof >= Uc.Size()) { + opserr << "PenaltyMP_FE::getResidual FATAL Error: Constrained DOF " << cdof << " out of bounds [0-" << Uc.Size() << "]\n"; + exit(-1); + } + UU(i) = Uc(cdof) - Uc0(i); + } + for (int i = 0; i < id2.Size(); ++i) { + int rdof = id2(i); + if (rdof < 0 || rdof >= Ur.Size()) { + opserr << "PenaltyMP_FE::getResidual FATAL Error: Retained DOF " << rdof << " out of bounds [0-" << Ur.Size() << "]\n"; + exit(-1); + } + UU(i+id1.Size()) = Ur(rdof) - Ur0(i); + } + + // compute residual + const Matrix& KK = getTangent(theNewIntegrator); + resid->addMatrixVector(0.0, KK, UU, -1.0); + + // done return *resid; } diff --git a/SRC/analysis/fe_ele/penalty/PenaltySP_FE.cpp b/SRC/analysis/fe_ele/penalty/PenaltySP_FE.cpp index 4cfe8cc9e6..dbd1ac9098 100644 --- a/SRC/analysis/fe_ele/penalty/PenaltySP_FE.cpp +++ b/SRC/analysis/fe_ele/penalty/PenaltySP_FE.cpp @@ -122,6 +122,7 @@ const Vector & PenaltySP_FE::getResidual(Integrator *theNewIntegrator) { double constraint = theSP->getValue(); + double initialValue = theSP->getInitialValue(); int constrainedDOF = theSP->getDOF_Number(); const Vector &nodeDisp = theNode->getTrialDisp(); @@ -135,7 +136,7 @@ PenaltySP_FE::getResidual(Integrator *theNewIntegrator) // is replace with the following to remove possible problems with // subtracting very small numbers - resid(0) = alpha * (constraint - nodeDisp(constrainedDOF)); + resid(0) = alpha * (constraint - (nodeDisp(constrainedDOF) - initialValue)); return resid; } diff --git a/SRC/analysis/handler/AutoConstraintHandler.cpp b/SRC/analysis/handler/AutoConstraintHandler.cpp new file mode 100644 index 0000000000..14bf7cdca7 --- /dev/null +++ b/SRC/analysis/handler/AutoConstraintHandler.cpp @@ -0,0 +1,532 @@ +/* ****************************************************************** ** +** 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.15 $ +// $Date: 2008-11-19 23:42:25 $ +// $Source: /usr/local/cvs/OpenSees/SRC/analysis/handler/AutoConstraintHandler.cpp,v $ + + +// Written: Massimo Petracca +// Created: June 2024. +// Revision: A +// +// What: "@(#) AutoConstraintHandler.C, revA" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +void* OPS_AutoConstraintHandler() +{ + // default parameters + bool verbose = false; + bool auto_penalty = true; + double auto_penalty_oom = 3.0; + double user_penalty = 0.0; + + // utils + const char* header = "constraints Auto <-verbose> <-autoPenalty $oom> <-userPenalty $userPenalty>"; + + // parse + bool auto_penalty_done = false; + bool user_penalty_done = false; + while (OPS_GetNumRemainingInputArgs() > 0) { + const char* type = OPS_GetString(); + if ((strcmp(type, "-verbose") == 0) || (strcmp(type, "-Verbose") == 0)) { + verbose = true; + } + else if ((strcmp(type, "-autoPenalty") == 0) || (strcmp(type, "-AutoPenalty") == 0)) { + if (OPS_GetNumRemainingInputArgs() == 0) { + opserr << "Error in: " << header << "\n"; + opserr << "$scale parameter not provided with -autoPenalty option\n"; + return 0; + } + int numData = 1; + if (OPS_GetDoubleInput(&numData, &auto_penalty_oom) != 0) { + opserr << "Error in: " << header << "\n"; + opserr << "Cannot get $oom parameter with -autoPenalty option\n"; + return 0; + } + if (user_penalty_done) { + opserr << "Error in: " << header << "\n"; + opserr << "Options -autoPenalty and -userPenalty are mutually exclusive\n"; + return 0; + } + auto_penalty = true; + auto_penalty_done = true; + } + else if ((strcmp(type, "-userPenalty") == 0) || (strcmp(type, "-UserPenalty") == 0)) { + if (OPS_GetNumRemainingInputArgs() == 0) { + opserr << "Error in: " << header << "\n"; + opserr << "$userPenalty parameter not provided with -userPenalty option\n"; + return 0; + } + int numData = 1; + if (OPS_GetDoubleInput(&numData, &user_penalty) != 0) { + opserr << "Error in: " << header << "\n"; + opserr << "Cannot get $userPenalty parameter with -userPenalty option\n"; + return 0; + } + if (auto_penalty_done) { + opserr << "Error in: " << header << "\n"; + opserr << "Options -autoPenalty and -userPenalty are mutually exclusive\n"; + return 0; + } + auto_penalty = false; + user_penalty_done = true; + } + } + + // done + return new AutoConstraintHandler( + verbose, + auto_penalty, + auto_penalty_oom, + user_penalty); +} + +namespace { + + /** + A utility class to compute, for each MP_Constraint, the penalty value according + to the (approximate) max stiffness value at the retained and constrained nodes. + */ + class PenaltyEvaluator + { + public: + PenaltyEvaluator() = delete; + /** + Full constructor. + + @param domain The domain + @param mps The list of MP_Constraints used for the auto-computation of penalty values + @param penalty_oom The order-of-magnitude to be added to the stiffness order-of-magnitude to obtain the penalty value + */ + PenaltyEvaluator(Domain* domain, const std::vector mps, double penalty_oom) + { + // store nodes involved in the input mps and init penalty to 0 + for (auto* mp : mps) { + if (mp) { + m_node_penalty[mp->getNodeRetained()] = 0.0; + m_node_penalty[mp->getNodeConstrained()] = 0.0; + } + } + // eval elements connected to those nodes + Element* elePtr; + ElementIter& theEles = domain->getElements(); + while ((elePtr = theEles()) != 0) { + if (!elePtr->isSubdomain()) { + // element nodes. process only if at least one node is involved + const ID& elenodes = elePtr->getExternalNodes(); + bool found = false; + for (int i = 0; i < elenodes.Size(); ++i) { + if (m_node_penalty.count(elenodes(i)) > 0) { + found = true; + break; + } + } + if (!found) continue; + // eval max diagonal entry for this element + const Matrix& K = elePtr->getInitialStiff(); + double kmax = 0.0; + for (int i = 0; i < K.noRows(); ++i) { + double ki = std::abs(K(i, i)); + kmax = std::max(kmax, ki); + } + // and accumulte the same to each node (it's just an approximation) + for (int i = 0; i < elenodes.Size(); ++i) { + auto iter = m_node_penalty.find(elenodes(i)); + if (iter != m_node_penalty.end()) + iter->second += kmax; + } + } + } + // now m_node_penalty contains the accumulated max diagonal stiffness entry + // at each node from each element connected to that node. + // now we can convert the stiffness value to the penalty value + for (auto& item : m_node_penalty) { + // the accumulated stiffness at this node + double kmax = item.second; + // its order of magnitude + double koom = std::round(std::log10(kmax)); + // compute the penalty value + item.second = std::pow(10.0, koom + penalty_oom); + } + } + /** + returns the penalty value for the input mp constraint + */ + inline double getPenaltyValue(const MP_Constraint* mp) const + { + double value = 0.0; + if (mp == nullptr) return value; + auto itr = m_node_penalty.find(mp->getNodeRetained()); + if (itr != m_node_penalty.end()) value = std::max(value, itr->second); + auto itc = m_node_penalty.find(mp->getNodeConstrained()); + if (itc != m_node_penalty.end()) value = std::max(value, itc->second); + return value; + } + + public: + // maps a node id to the penalty value at that node + std::unordered_map m_node_penalty; + }; + +} + +AutoConstraintHandler::AutoConstraintHandler() + : ConstraintHandler(HANDLER_TAG_AutoConstraintHandler) +{ + +} + +AutoConstraintHandler::AutoConstraintHandler( + bool _verbose, + bool _auto_penalty, + double _auto_penalty_oom, + double _user_penalty) + : ConstraintHandler(HANDLER_TAG_AutoConstraintHandler) + , verbose(_verbose) + , auto_penalty(_auto_penalty) + , auto_penalty_oom(_auto_penalty_oom) + , user_penalty(_user_penalty) +{ +} + +AutoConstraintHandler::~AutoConstraintHandler() +{ + +} + +int +AutoConstraintHandler::handle(const ID* nodesLast) +{ + // first check links exist to a Domain and an AnalysisModel object + Domain* theDomain = this->getDomainPtr(); + AnalysisModel* theModel = this->getAnalysisModelPtr(); + Integrator* theIntegrator = this->getIntegratorPtr(); + if ((theDomain == 0) || (theModel == 0) || (theIntegrator == 0)) { + opserr << "WARNING AutoConstraintHandler::handle() - "; + opserr << " setLinks() has not been called\n"; + return -1; + } + + // create a multimap (key = NodeID, value = SP_Constraint). + // multimap allows for duplicate keys because we may have multiple SP + // constraints on the same node. + // in this way we have : + // a) a fast access to the key (to tell whether a node is constrained or not) + // b) and a range of SP_Constraints acting on that node + std::unordered_multimap sp_map; + { + SP_ConstraintIter& theSPs = theDomain->getDomainAndLoadPatternSPs(); + SP_Constraint* theSP; + while ((theSP = theSPs()) != 0) + sp_map.emplace(std::make_pair(theSP->getNodeTag(), theSP)); + } + + // create a DOF_Group for each Node and add it to the AnalysisModel + // create a TrasformationDOF_Group for nodes constrained in SP_Constraints + int countDOF = 0; + { + NodeIter& theNodes = theDomain->getNodes(); + Node* nodPtr; + int numDofGrp = 0; + while ((nodPtr = theNodes()) != 0) { + + // process this node + DOF_Group* dofPtr = 0; + int nodeTag = nodPtr->getTag(); + + // if the node is constrained in an SP constraint + // handle it with the transformation method + auto it_range = sp_map.equal_range(nodeTag); + if (it_range.first != it_range.second) { + + // this node is constrained in 1 or more SP_Constraints. + // handle it with a TransformationDOF_Group. + // note: the TransformationConstraintHandler pointer is useless in TDOF... + TransformationDOF_Group* tDofPtr = + new TransformationDOF_Group(numDofGrp++, nodPtr, nullptr); + dofPtr = tDofPtr; + + // count and add all SP_Constraints. + // if verbose, keep track of conflicting constraints + int numSPs = 0; + for (auto it = it_range.first; it != it_range.second; it++) { + tDofPtr->addSP_Constraint(*(it->second)); + numSPs++; + } + + // add the DOF to the array + theDOFs.push_back(tDofPtr); + // count only free DOFs + countDOF += nodPtr->getNumberDOF() - numSPs; + + } + else { + + // free node, create an ordinary DOF_Group + dofPtr = new DOF_Group(numDofGrp++, nodPtr); + // count all DOFs + countDOF += nodPtr->getNumberDOF(); + } + + // add the DOF Group to the model + nodPtr->setDOF_GroupPtr(dofPtr); + theModel->addDOF_Group(dofPtr); + } + } + + // set the number of equations + theModel->setNumEqn(countDOF); + + // now see if we have to set any of the dof's to -3 + int count3 = 0; + if (nodesLast != 0) { + for (int i = 0; i < nodesLast->Size(); i++) { + int nodeID = (*nodesLast)(i); + Node* nodPtr = theDomain->getNode(nodeID); + if (nodPtr != 0) { + DOF_Group* dofPtr = nodPtr->getDOF_GroupPtr(); + const ID& id = dofPtr->getID(); + // set all the dof values to -3 + for (int j = 0; j < id.Size(); j++) { + if (id(j) == -2) { + dofPtr->setID(j, -3); + count3++; + } + else { + opserr << "WARNING AutoConstraintHandler::handle()" + " - boundary sp constraint in subdomain this should not be " + "- results suspect \n"; + } + } + } + } + } + + // create the FE_Elements for the Elements and add to the AnalysisModel + ElementIter& theEle = theDomain->getElements(); + Element* elePtr; + int numFeEle = 0; + FE_Element* fePtr; + + // first standard elements + while ((elePtr = theEle()) != 0) { + // only create an FE_Element for a subdomain element if it does not + // do independent analysis .. then subdomain part of this analysis so create + // an FE_element & set subdomain to point to it. + if (elePtr->isSubdomain() == true) { + Subdomain* theSub = (Subdomain*)elePtr; + if (theSub->doesIndependentAnalysis() == false) { + fePtr = new FE_Element(numFeEle++, elePtr); + theModel->addFE_Element(fePtr); + theSub->setFE_ElementPtr(fePtr); + } + } + else { + // just a regular element .. create an FE_Element for it & add to AnalysisModel + fePtr = new FE_Element(numFeEle++, elePtr); + theModel->addFE_Element(fePtr); + } + } + + // then MP Penalty elements + // Note: now we get all MPs and handle them with penalty. + // In the future we may split valid MPs and conflicting MPs. valid one will be handled + // by the transformation method + std::vector mps_penalty; + { + MP_ConstraintIter& theMPs = theDomain->getMPs(); + MP_Constraint* mpPtr; + while ((mpPtr = theMPs()) != 0) + mps_penalty.push_back(mpPtr); + + } + std::shared_ptr peval; + if (auto_penalty) + peval = std::make_shared(theDomain, mps_penalty, auto_penalty_oom); + for(MP_Constraint* mp : mps_penalty) { + double penalty = auto_penalty ? peval->getPenaltyValue(mp) : user_penalty; + fePtr = new PenaltyMP_FE(numFeEle, *theDomain, *mp, penalty); + theModel->addFE_Element(fePtr); + numFeEle++; + } + + // give some info + if (verbose) { + std::stringstream ss; + ss << "AutoConstraint Handler report:\n"; + ss << "+ SP Constraints:\n"; + ss << " + " << sp_map.size() << " constraints handled with the Transformation method\n"; + if (sp_map.size() > 0) { + NodeIter& theNodes = theDomain->getNodes(); + Node* nodPtr; + int numDofGrp = 0; + while ((nodPtr = theNodes()) != 0) { + auto it_range = sp_map.equal_range(nodPtr->getTag()); + if (it_range.first == it_range.second) continue; + std::vector sp_per_dof(static_cast(nodPtr->getNumberDOF()), 0); + for (auto it = it_range.first; it != it_range.second; it++) { + SP_Constraint* sp = it->second; + sp_per_dof[static_cast(sp->getDOF_Number())] += 1; + } + ss << " - Node: " << nodPtr->getTag() << " [ "; + for (auto i : sp_per_dof) ss << i << " "; + ss << "]\n"; + for (auto it = it_range.first; it != it_range.second; it++) { + SP_Constraint* sp = it->second; + int how_many = sp_per_dof[static_cast(sp->getDOF_Number())]; + if (how_many > 1) { + ss << " ! Warning: @DOF " + << sp->getDOF_Number() + << ". Found " << how_many + << " SP constraints. Only the SP (tag = " + << sp->getTag() << ", value = " << sp->getValue() + << ") will be imposed\n"; + } + } + } + } + ss << "+ MP Constraints:\n"; + ss << " + " << mps_penalty.size() << " constraints handled with the Penalty method\n"; + if (auto_penalty) { + ss << " + Automatic Penalty values for each MP Constraint:\n"; + for (MP_Constraint* mp : mps_penalty) + ss << " + MP(tag = " << mp->getTag() << ") = " << std::scientific << peval->getPenaltyValue(mp) << "\n"; + } + else { + ss << " + Uniform User-Defined penalty = " << std::scientific << user_penalty << "\n"; + } + std::string ss_value = ss.str(); + opserr << ss_value.c_str(); + } + + // done + return count3; +} + +void +AutoConstraintHandler::clearAll(void) +{ + // reset the TransformationDOF_Group vector + // don't delete the pointer inside the vector (owned by the model) + theDOFs.clear(); + + // for the nodes reset the DOF_Group pointers to 0 + Domain* theDomain = this->getDomainPtr(); + if (theDomain == 0) + return; + NodeIter& theNod = theDomain->getNodes(); + Node* nodPtr; + while ((nodPtr = theNod()) != 0) + nodPtr->setDOF_GroupPtr(0); +} + +int +AutoConstraintHandler::sendSelf(int cTag, Channel& theChannel) +{ + Vector data(4); + int result = 0; + data(0) = static_cast(verbose); + data(1) = static_cast(auto_penalty); + data(2) = auto_penalty_oom; + data(3) = user_penalty; + result = theChannel.sendVector(this->getDbTag(), cTag, data); + if (result != 0) + opserr << "AutoConstraintHandler::sendSelf() - error sending Vector\n"; + return result; +} + +int +AutoConstraintHandler::recvSelf(int cTag, + Channel& theChannel, + FEM_ObjectBroker& theBroker) +{ + Vector data(4); + int result = 0; + result = theChannel.recvVector(this->getDbTag(), cTag, data); + verbose = static_cast(data(0)); + auto_penalty = static_cast(data(1)); + auto_penalty_oom = data(2); + user_penalty = data(3); + if (result != 0) + opserr << "AutoConstraintHandler::recvSelf() - error receiving Vector\n"; + return result; +} + +int +AutoConstraintHandler::applyLoad(void) +{ + // enforse SP constraints + // is there a reason for a backward loop? + std::size_t numDOF = theDOFs.size(); + for (std::size_t i = 1; i <= numDOF; i++) + theDOFs[numDOF - i]->enforceSPs(1); + for (std::size_t i = 1; i <= numDOF; i++) + theDOFs[numDOF - i]->enforceSPs(0); + return 0; +} + +int +AutoConstraintHandler::doneNumberingDOF(void) +{ + // iterate through the DOF_Groups telling them that their ID has now been set + DOF_GrpIter& theDOFS = this->getAnalysisModelPtr()->getDOFs(); + DOF_Group* dofPtr; + while ((dofPtr = theDOFS()) != 0) + dofPtr->doneID(); + + // iterate through the FE_Element getting them to set their IDs + // done using base class implementation + return ConstraintHandler::doneNumberingDOF(); +} diff --git a/SRC/analysis/handler/AutoConstraintHandler.h b/SRC/analysis/handler/AutoConstraintHandler.h new file mode 100644 index 0000000000..f2de41f613 --- /dev/null +++ b/SRC/analysis/handler/AutoConstraintHandler.h @@ -0,0 +1,88 @@ +/* ****************************************************************** ** +** 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.4 $ +// $Date: 2005-11-28 21:35:54 $ +// $Source: /usr/local/cvs/OpenSees/SRC/analysis/handler/AutoConstraintHandler.h,v $ + + +// Written: Massimo Petracca +// Created: June 2024. +// Revision: A +// +// Description: This file contains the class definition for +// AutoConstraintHandler. AutoConstraintHandler is a +// constraint handler for handling constraints using the different methods. +// +// For each element and degree-of-freedom at a node it constructs: +// +// 1) regular FE_Element and DOF_Groups if there is no SP_Constraint or MP_Constraint; +// 2) TransformationDOF_Group for SP constraints (as in the Transformation method) +// 3) PenaltyMP_FE for MP constraints +// +// Notes: +// 1) For each PenaltyMP_FE, by default it automatically selects a proper penalty +// value based on the stiffness values found on the DOFs involved in the +// MP constraint +// +// What: "@(#) AutoConstraintHandler.h, revA" + +#ifndef AutoConstraintHandler_h +#define AutoConstraintHandler_h + +#include +#include + +class FE_Element; +class TransformationDOF_Group; + +class AutoConstraintHandler : public ConstraintHandler +{ +public: + AutoConstraintHandler(); + AutoConstraintHandler( + bool _verbose, + bool _auto_penalty, + double _auto_penalty_oom, + double _user_penalty); + ~AutoConstraintHandler(); + + int handle(const ID* nodesNumberedLast = 0); + int applyLoad(); + void clearAll(void); + int doneNumberingDOF(void); + + int sendSelf(int commitTag, Channel& theChannel); + int recvSelf(int commitTag, Channel& theChannel, + FEM_ObjectBroker& theBroker); + +private: + bool verbose = false; + bool auto_penalty = true; + double auto_penalty_oom = 3.0; + double user_penalty = 0.0; + std::vector theDOFs; +}; + +#endif + + + + diff --git a/SRC/analysis/handler/CMakeLists.txt b/SRC/analysis/handler/CMakeLists.txt index 7053dd1de2..1bf8a5da14 100644 --- a/SRC/analysis/handler/CMakeLists.txt +++ b/SRC/analysis/handler/CMakeLists.txt @@ -11,12 +11,14 @@ target_sources(OPS_Analysis PenaltyConstraintHandler.cpp LagrangeConstraintHandler.cpp TransformationConstraintHandler.cpp + AutoConstraintHandler.cpp PUBLIC ConstraintHandler.h PlainHandler.h PenaltyConstraintHandler.h LagrangeConstraintHandler.h TransformationConstraintHandler.h + AutoConstraintHandler.h ) target_include_directories(OPS_Analysis PUBLIC ${CMAKE_CURRENT_LIST_DIR}) diff --git a/SRC/analysis/handler/Makefile b/SRC/analysis/handler/Makefile index de8e94a0e0..c2144ab114 100644 --- a/SRC/analysis/handler/Makefile +++ b/SRC/analysis/handler/Makefile @@ -2,6 +2,7 @@ include ../../../Makefile.def OBJS = ConstraintHandler.o PlainHandler.o \ PenaltyConstraintHandler.o LagrangeConstraintHandler.o \ + AutoConstraintHandler.o \ TransformationConstraintHandler.o PenaltyHandlerNoHomoSPMultipliers.o # Compilation control diff --git a/SRC/classTags.h b/SRC/classTags.h index 2527b6910e..1886497d76 100644 --- a/SRC/classTags.h +++ b/SRC/classTags.h @@ -977,6 +977,7 @@ #define HANDLER_TAG_PenaltyConstraintHandler 3 #define HANDLER_TAG_TransformationConstraintHandler 4 #define HANDLER_TAG_PenaltyHandlerNoHomoSPMultipliers 5 +#define HANDLER_TAG_AutoConstraintHandler 6 #define NUMBERER_TAG_DOF_Numberer 1 #define NUMBERER_TAG_PlainNumberer 2 diff --git a/SRC/domain/constraints/MP_Constraint.cpp b/SRC/domain/constraints/MP_Constraint.cpp index 550a5578e5..93447645d3 100644 --- a/SRC/domain/constraints/MP_Constraint.cpp +++ b/SRC/domain/constraints/MP_Constraint.cpp @@ -41,6 +41,7 @@ #include #include #include +#include static int numMPs = 0; static int nextTag = 0; @@ -222,6 +223,12 @@ MP_Constraint::MP_Constraint(int nodeRetain, int nodeConstr, opserr << "MP_Constraint::MP_Constraint - ran out of memory 1\n"; exit(-1); } + + // resize initial state + Uc0.resize(constrDOF->Size()); + Uc0.Zero(); + Ur0.resize(retainDOF->Size()); + Ur0.Zero(); } @@ -246,6 +253,12 @@ MP_Constraint::MP_Constraint(int nodeRetain, int nodeConstr, Matrix &constr, opserr << "MP_Constraint::MP_Constraint - ran out of memory 2\n"; exit(-1); } + + // resize initial state + Uc0.resize(constrDOF->Size()); + Uc0.Zero(); + Ur0.resize(retainDOF->Size()); + Ur0.Zero(); } @@ -263,9 +276,46 @@ MP_Constraint::~MP_Constraint() numMPs--; if (numMPs == 0) nextTag = 0; +} + +void MP_Constraint::setDomain(Domain* theDomain) +{ + // store initial state + if (theDomain) { + Node* theRetainedNode = theDomain->getNode(nodeRetained); + Node* theConstrainedNode = theDomain->getNode(nodeConstrained); + if (theRetainedNode == 0 || theConstrainedNode == 0) { + opserr << "FATAL MP_Constraint::setDomain() - Constrained or Retained"; + opserr << " Node does not exist in Domain\n"; + opserr << nodeRetained << " " << nodeConstrained << endln; + exit(-1); + } + const Vector& Uc = theConstrainedNode->getTrialDisp(); + const Vector& Ur = theRetainedNode->getTrialDisp(); + const ID& idc = getConstrainedDOFs(); + const ID& idr = getRetainedDOFs(); + for (int i = 0; i < idc.Size(); ++i) { + int cdof = idc(i); + if (cdof < 0 || cdof >= Uc.Size()) { + opserr << "MP_Constraint::setDomain FATAL Error: Constrained DOF " << cdof << " out of bounds [0-" << Uc.Size() << "]\n"; + exit(-1); + } + Uc0(i) = Uc(cdof); + } + for (int i = 0; i < idr.Size(); ++i) { + int rdof = idr(i); + if (rdof < 0 || rdof >= Ur.Size()) { + opserr << "MP_Constraint::setDomain FATAL Error: Retained DOF " << rdof << " out of bounds [0-" << Ur.Size() << "]\n"; + exit(-1); + } + Ur0(i) = Ur(rdof); + } + } + + // call base class implementation + DomainComponent::setDomain(theDomain); } - int MP_Constraint::getNodeRetained(void) const { @@ -332,6 +382,16 @@ MP_Constraint::getConstraint(void) // return the constraint matrix Ccr return *constraint; +} + +const Vector& MP_Constraint::getConstrainedDOFsInitialDisplacement(void) const +{ + return Uc0; +} + +const Vector& MP_Constraint::getRetainedDOFsInitialDisplacement(void) const +{ + return Ur0; } int @@ -391,6 +451,26 @@ MP_Constraint::sendSelf(int cTag, Channel &theChannel) } } + // send initial displacement vectors. + // we need 2 database tags because they have the same size, + // but we can reuse the tags used for ID objects, since they go into different files + if (Uc0.Size() > 0) { + int result = theChannel.sendVector(dbTag1, cTag, Uc0); + if (result < 0) { + opserr << "WARNING MP_Constraint::sendSelf "; + opserr << "- error sending constrained initial displacement\n"; + return result; + } + } + if (Ur0.Size() > 0) { + int result = theChannel.sendVector(dbTag2, cTag, Ur0); + if (result < 0) { + opserr << "WARNING MP_Constraint::sendSelf "; + opserr << "- error sending retained initial displacement\n"; + return result; + } + } + return 0; } @@ -442,12 +522,40 @@ MP_Constraint::recvSelf(int cTag, Channel &theChannel, retainDOF = new ID(size); int result = theChannel.recvID(dbTag2, cTag, *retainDOF); if (result < 0) { - opserr << "WARNING MP_Retainaint::recvSelf "; + opserr << "WARNING MP_Constraint::recvSelf "; opserr << "- error receiving retained data\n"; return result; } } + // recv initial displacement vectors. + // we need 2 database tags because they have the same size, + // but we can reuse the tags used for ID objects, since they go into different files + if (constrDOF && constrDOF->Size() > 0) + Uc0.resize(constrDOF->Size()); + else + Uc0 = Vector(); + if (retainDOF && retainDOF->Size() > 0) + Ur0.resize(retainDOF->Size()); + else + Ur0 = Vector(); + if (Uc0.Size() > 0) { + int result = theChannel.recvVector(dbTag1, cTag, Uc0); + if (result < 0) { + opserr << "WARNING MP_Constraint::recvSelf "; + opserr << "- error receiving constrained initial displacement\n"; + return result; + } + } + if (Ur0.Size() > 0) { + int result = theChannel.recvVector(dbTag2, cTag, Ur0); + if (result < 0) { + opserr << "WARNING MP_Constraint::recvSelf "; + opserr << "- error receiving retained initial displacement\n"; + return result; + } + } + return 0; } diff --git a/SRC/domain/constraints/MP_Constraint.h b/SRC/domain/constraints/MP_Constraint.h index 417405b1cc..83ce3ed3b8 100644 --- a/SRC/domain/constraints/MP_Constraint.h +++ b/SRC/domain/constraints/MP_Constraint.h @@ -68,6 +68,9 @@ class MP_Constraint : public DomainComponent // destructor virtual ~MP_Constraint(); + // domain component + void setDomain(Domain* theDomain); + // method to get information about the constraint virtual int getNodeRetained(void) const; virtual int getNodeConstrained(void) const; @@ -76,6 +79,8 @@ class MP_Constraint : public DomainComponent virtual int applyConstraint(double pseudoTime); virtual bool isTimeVarying(void) const; virtual const Matrix &getConstraint(void); + virtual const Vector &getConstrainedDOFsInitialDisplacement(void) const; + virtual const Vector &getRetainedDOFsInitialDisplacement(void) const; // methods for output virtual int sendSelf(int commitTag, Channel &theChannel); @@ -92,7 +97,8 @@ class MP_Constraint : public DomainComponent Matrix *constraint; // pointer to the constraint matrix ID *constrDOF; // ID of constrained DOF at constrained node ID *retainDOF; // ID of related DOF at retained node - + Vector Uc0; // initial displacement at constrained DOFs (same size as constrDOF) + Vector Ur0; // initial displacement at retained node (same size as retainDOF) int dbTag1, dbTag2; // need a dbTag for the two ID's }; diff --git a/SRC/domain/constraints/SP_Constraint.cpp b/SRC/domain/constraints/SP_Constraint.cpp index 8a79b96ddc..0c30e18a08 100644 --- a/SRC/domain/constraints/SP_Constraint.cpp +++ b/SRC/domain/constraints/SP_Constraint.cpp @@ -257,7 +257,7 @@ int SP_Constraint_GetNextTag(void) { // constructor for FEM_ObjectBroker SP_Constraint::SP_Constraint(int clasTag) :DomainComponent(0,clasTag), - nodeTag(0), dofNumber(0), valueR(0.0), valueC(0.0), isConstant(true), + nodeTag(0), dofNumber(0), valueR(0.0), valueC(0.0), initialValue(0.0), isConstant(true), loadPatternTag(-1) { numSPs++; @@ -266,7 +266,7 @@ SP_Constraint::SP_Constraint(int clasTag) // constructor for a subclass to use SP_Constraint::SP_Constraint(int node, int ndof, int clasTag) :DomainComponent(nextTag++, clasTag), - nodeTag(node), dofNumber(ndof), valueR(0.0), valueC(0.0), isConstant(true), + nodeTag(node), dofNumber(ndof), valueR(0.0), valueC(0.0), initialValue(0.0), isConstant(true), loadPatternTag(-1) // valueC is set to 1.0 so that homo will be false when recvSelf() invoked // should be ok as valueC cannot be used by subclasses and subclasses should @@ -278,7 +278,7 @@ SP_Constraint::SP_Constraint(int node, int ndof, int clasTag) // constructor for object of type SP_Constraint SP_Constraint::SP_Constraint(int node, int ndof, double value, bool ISconstant) :DomainComponent(nextTag++, CNSTRNT_TAG_SP_Constraint), - nodeTag(node), dofNumber(ndof), valueR(value), valueC(value), isConstant(ISconstant), + nodeTag(node), dofNumber(ndof), valueR(value), valueC(value), initialValue(0.0), isConstant(ISconstant), loadPatternTag(-1) { numSPs++; @@ -313,6 +313,13 @@ SP_Constraint::getValue(void) return valueC; } +double +SP_Constraint::getInitialValue(void) +{ + // return the initial value of the constraint + return initialValue; +} + int SP_Constraint::applyConstraint(double loadFactor) { @@ -345,10 +352,34 @@ SP_Constraint::getLoadPatternTag(void) const return loadPatternTag; } +void +SP_Constraint::setDomain(Domain* theDomain) +{ + // store initial state + if (theDomain) { + Node* theNode = theDomain->getNode(nodeTag); + if (theNode == 0) { + opserr << "FATAL SP_Constraint::setDomain() - Constrained"; + opserr << " Node does not exist in Domain\n"; + opserr << nodeTag << endln; + exit(-1); + } + const Vector& U = theNode->getTrialDisp(); + if (dofNumber < 0 || dofNumber >= U.Size()) { + opserr << "SP_Constraint::setDomain FATAL Error: Constrained DOF " << dofNumber << " out of bounds [0-" << U.Size() << "]\n"; + exit(-1); + } + initialValue = U(dofNumber); + } + + // call base class implementation + DomainComponent::setDomain(theDomain); +} + int SP_Constraint::sendSelf(int cTag, Channel &theChannel) { - static Vector data(8); // we send as double to avoid having + static Vector data(9); // we send as double to avoid having // to send two messages. data(0) = this->getTag(); data(1) = nodeTag; @@ -362,6 +393,7 @@ SP_Constraint::sendSelf(int cTag, Channel &theChannel) data(6) = this->getLoadPatternTag(); data(7) = nextTag; + data(8) = initialValue; int result = theChannel.sendVector(this->getDbTag(), cTag, data); if (result != 0) { @@ -376,7 +408,7 @@ int SP_Constraint::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker) { - static Vector data(8); // we sent the data as double to avoid having to send + static Vector data(9); // we sent the data as double to avoid having to send // two messages int result = theChannel.recvVector(this->getDbTag(), cTag, data); if (result < 0) { @@ -399,6 +431,7 @@ SP_Constraint::recvSelf(int cTag, Channel &theChannel, this->setLoadPatternTag((int)data(6)); nextTag = (int)data(7); + initialValue = data(8); return 0; } @@ -409,7 +442,7 @@ SP_Constraint::Print(OPS_Stream &s, int flag) { s << "SP_Constraint: " << this->getTag(); s << "\t Node: " << nodeTag << " DOF: " << dofNumber+1; - s << " ref value: " << valueR << " current value: " << valueC << endln; + s << " ref value: " << valueR << " current value: " << valueC << " initial value: " << initialValue << endln; } diff --git a/SRC/domain/constraints/SP_Constraint.h b/SRC/domain/constraints/SP_Constraint.h index c36945d092..cd2c656bb3 100644 --- a/SRC/domain/constraints/SP_Constraint.h +++ b/SRC/domain/constraints/SP_Constraint.h @@ -56,10 +56,13 @@ class SP_Constraint : public DomainComponent virtual int getDOF_Number(void) const; virtual int applyConstraint(double loadFactor); virtual double getValue(void); + virtual double getInitialValue(void); virtual bool isHomogeneous(void) const; virtual void setLoadPatternTag(int loadPaternTag); virtual int getLoadPatternTag(void) const; + void setDomain(Domain* theDomain); + virtual int sendSelf(int commitTag, Channel &theChannel); virtual int recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker); @@ -72,6 +75,7 @@ class SP_Constraint : public DomainComponent double valueR; // the reference value double valueC; // if constant = the reference value, if not constant = // the reference value * load factor + double initialValue; // the value of the dof when the sp constrain is added to the domain bool isConstant; // flag indicating if constant int loadPatternTag; }; diff --git a/SRC/interpreter/OpenSeesCommands.cpp b/SRC/interpreter/OpenSeesCommands.cpp index d8358f3f9f..c5b7fd62c0 100644 --- a/SRC/interpreter/OpenSeesCommands.cpp +++ b/SRC/interpreter/OpenSeesCommands.cpp @@ -1531,6 +1531,9 @@ int OPS_ConstraintHandler() } else if (strcmp(type,"Transformation") == 0) { theHandler = (ConstraintHandler*)OPS_TransformationConstraintHandler(); + } else if (strcmp(type,"Auto") == 0) { + theHandler = (ConstraintHandler*)OPS_AutoConstraintHandler(); + } else { opserr<<"WARNING unknown ConstraintHandler type "< #include +extern void* OPS_AutoConstraintHandler(void); + // numberers #include #include @@ -3779,6 +3781,13 @@ specifyConstraintHandler(ClientData clientData, Tcl_Interp *interp, int argc, theHandler = new TransformationConstraintHandler(); } + else if (strcmp(argv[1],"Auto") == 0) { + OPS_ResetInputNoBuilder(clientData, interp, 2, argc, argv, &theDomain); + theHandler = (ConstraintHandler*)OPS_AutoConstraintHandler(); + if (theHandler == 0) + return TCL_ERROR; + } + else { opserr << "WARNING No ConstraintHandler type exists (Plain, Penalty,\n"; opserr << " Lagrange, Transformation) only\n"; diff --git a/Win32/proj/analysis/analysis.vcxproj b/Win32/proj/analysis/analysis.vcxproj index 73d0a5a708..228ce29378 100644 --- a/Win32/proj/analysis/analysis.vcxproj +++ b/Win32/proj/analysis/analysis.vcxproj @@ -214,6 +214,7 @@ + @@ -329,6 +330,7 @@ + diff --git a/Win32/proj/analysis/analysis.vcxproj.filters b/Win32/proj/analysis/analysis.vcxproj.filters index 21bed9aa3c..4652c1361a 100644 --- a/Win32/proj/analysis/analysis.vcxproj.filters +++ b/Win32/proj/analysis/analysis.vcxproj.filters @@ -320,6 +320,9 @@ handler + + handler + lineSearch @@ -661,6 +664,9 @@ handler + + handler + lineSearch diff --git a/Win64/proj/analysis/analysis.vcxproj b/Win64/proj/analysis/analysis.vcxproj index 410da1ce1b..ec79f0733e 100644 --- a/Win64/proj/analysis/analysis.vcxproj +++ b/Win64/proj/analysis/analysis.vcxproj @@ -306,6 +306,7 @@ + @@ -422,6 +423,7 @@ + diff --git a/Win64/proj/analysis/analysis.vcxproj.filters b/Win64/proj/analysis/analysis.vcxproj.filters index de96a88243..a477859334 100644 --- a/Win64/proj/analysis/analysis.vcxproj.filters +++ b/Win64/proj/analysis/analysis.vcxproj.filters @@ -320,6 +320,9 @@ handler + + handler + lineSearch @@ -670,6 +673,9 @@ handler + + handler + lineSearch