Merge pull request OpenSees#1450 from MassimoPetracca/constraint-hand…

…ler-incremental-form

Improving constraint handlers

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

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 \

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 ";

SRC/analysis/dof_grp/LagrangeDOF_Group.cpp

@@ -42,14 +42,23 @@
 #include <SP_Constraint.h>
 #include <MP_Constraint.h>
 
+//#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  

SRC/analysis/dof_grp/LagrangeDOF_Group.h

@@ -39,6 +39,7 @@
 // What: "@(#) LagrangeDOF_Group.h, revA"
 
 #include <DOF_Group.h>
+#include <Vector.h>
 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;
 
 };
 

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; i<numDOF; i++) {
     if (theSPs[i] != 0)
+#ifdef TRANSF_INCREMENTAL_MP
+      (*unbalance)(i) = 0.0; // don't enfore the SP here as in incrNodeDisp!
+#else
       (*unbalance)(i) = disp(i);
+#endif // TRANSF_INCREMENTAL_MP
   }
+
+#ifdef TRANSF_INCREMENTAL_MP
+  myNode->incrTrialDisp(*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; i<numDof; i++)
       if (theSPs[i] != 0) {
 	double value = theSPs[i]->getValue();
-	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
       }
   } 
 

SRC/analysis/dof_grp/TransformationDOF_Group.h

@@ -40,6 +40,15 @@
 
 #include <DOF_Group.h>
 
+// 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 <Vector.h>
+#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

SRC/analysis/fe_ele/lagrange/LagrangeMP_FE.cpp

@@ -46,11 +46,17 @@
 #include <MP_Constraint.h>
 #include <DOF_Group.h>
 
+// 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;
 }