My title

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.1.1.1 $
00022 // $Date: 2000/09/15 08:23:16 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/fe_ele/penalty/PenaltyMP_FE.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/analysis/fe_ele/penalty/PenaltyMP_FE.C
00027 //
00028 // Written: fmk 
00029 // Created: 11/96
00030 // Revision: A
00031 //
00032 // Purpose: This file contains the code for implementing the methods
00033 // of the PenaltyMP_FE class interface.
00034 //
00035 // the interface:
00036 
00037 #include <PenaltyMP_FE.h>
00038 #include <stdlib.h>
00039 
00040 #include <Element.h>
00041 #include <Domain.h>
00042 #include <Node.h>
00043 #include <DOF_Group.h>
00044 #include <Integrator.h>
00045 #include <Subdomain.h>
00046 #include <AnalysisModel.h>
00047 #include <Matrix.h>
00048 #include <Vector.h>
00049 #include <Node.h>
00050 #include <MP_Constraint.h>
00051 #include <DOF_Group.h>
00052 
00053 PenaltyMP_FE::PenaltyMP_FE(Domain &theDomain, 
00054       MP_Constraint &TheMP, double Alpha)
00055 :FE_Element(2,(TheMP.getConstrainedDOFs()).Size()+
00056  (TheMP.getRetainedDOFs()).Size()),
00057  theMP(&TheMP), theConstrainedNode(0) , theRetainedNode(0),
00058  tang(0), resid(0), C(0), alpha(Alpha)
00059 {
00060     
00061     int size;
00062     const ID &id1 = theMP->getConstrainedDOFs();
00063     size = id1.Size();
00064     const ID &id2 = theMP->getRetainedDOFs();    
00065     size += id2.Size();
00066 
00067     tang = new Matrix(size,size);
00068     resid = new Vector(size);
00069     C = new Matrix(id1.Size(),size);
00070 
00071     if (tang == 0 || resid == 0 || C == 0 ||
00072  tang->noCols() != size || C->noCols() != size || 
00073  resid->Size() != size) {
00074  cerr << "FATAL PenaltyMP_FE::PenaltyMP_FE() - out of memory\n";
00075  exit(-1);
00076     }
00077      
00078     theRetainedNode = theDomain.getNode(theMP->getNodeRetained());    
00079     theConstrainedNode = theDomain.getNode(theMP->getNodeConstrained());
00080 
00081     if (theRetainedNode == 0 || theConstrainedNode == 0) {
00082  cerr << "FATAL PenaltyMP_FE::PenaltyMP_FE() - Constrained or Retained";
00083  cerr << " Node does not exist in Domain\n";
00084  exit(-1);
00085     } 
00086 
00087 
00088     // set up the dof groups tags
00089     DOF_Group *dofGrpPtr = 0;
00090     dofGrpPtr = theRetainedNode->getDOF_GroupPtr();
00091     if (dofGrpPtr != 0) 
00092  myDOF_Groups(0) = dofGrpPtr->getTag();     
00093     else 
00094  cerr << "WARNING PenaltyMP_FE::PenaltyMP_FE() - node no Group yet?\n"; 
00095     dofGrpPtr = theConstrainedNode->getDOF_GroupPtr();
00096     if (dofGrpPtr != 0) 
00097  myDOF_Groups(1) = dofGrpPtr->getTag();         
00098     else
00099  cerr << "WARNING PenaltyMP_FE::PenaltyMP_FE() - node no Group yet?\n"; 
00100     
00101     
00102     if (theMP->isTimeVarying() == false) {
00103  this->determineTangent();
00104  // we can free up the space taken by C as it is no longer needed
00105  if (C != 0)
00106      delete C;
00107  C = 0;
00108     }
00109 }
00110 
00111 PenaltyMP_FE::~PenaltyMP_FE()
00112 {
00113     if (tang != 0)
00114  delete tang;
00115     if (resid != 0)
00116  delete resid;
00117     if (C != 0)
00118  delete C;
00119 }    
00120 
00121 // void setID(int index, int value);
00122 // Method to set the correMPonding index of the ID to value.
00123 int
00124 PenaltyMP_FE::setID(void)
00125 {
00126     int result = 0;
00127 
00128     // first determine the IDs in myID for those DOFs marked
00129     // as constrained DOFs, this is obtained from the DOF_Group
00130     // associated with the constrained node
00131     DOF_Group *theConstrainedNodesDOFs = theConstrainedNode->getDOF_GroupPtr();
00132     if (theConstrainedNodesDOFs == 0) {
00133  cerr << "WARNING PenaltyMP_FE::setID(void)";
00134  cerr << " - no DOF_Group with Constrained Node\n";
00135  return -2;
00136     }    
00137 
00138     const ID &constrainedDOFs = theMP->getConstrainedDOFs();
00139     const ID &theConstrainedNodesID = theConstrainedNodesDOFs->getID();    
00140     
00141     int size1 = constrainedDOFs.Size();
00142     for (int i=0; i<size1; i++) {
00143  int constrained = constrainedDOFs(i);
00144  if (constrained < 0 || 
00145      constrained >= theConstrainedNode->getNumberDOF()) {
00146      
00147      cerr << "WARNING PenaltyMP_FE::setID(void) - unknown DOF ";
00148      cerr << constrained << " at Node\n";
00149      myID(i) = -1; // modify so nothing will be added to equations
00150      result = -3;
00151  }     
00152  else {
00153      if (constrained >= theConstrainedNodesID.Size()) {
00154   cerr << "WARNING PenaltyMP_FE::setID(void) - ";
00155   cerr << " Nodes DOF_Group too small\n";
00156   myID(i) = -1; // modify so nothing will be added to equations
00157   result = -4;
00158      }
00159      else
00160   myID(i) = theConstrainedNodesID(constrained);
00161  }
00162     }
00163     
00164     // now determine the IDs for the retained dof's
00165     DOF_Group *theRetainedNodesDOFs = theRetainedNode->getDOF_GroupPtr();
00166     if (theRetainedNodesDOFs == 0) {
00167  cerr << "WARNING PenaltyMP_FE::setID(void)";
00168  cerr << " - no DOF_Group with Retained Node\n";
00169  return -2;
00170     }    
00171     
00172     const ID &RetainedDOFs = theMP->getRetainedDOFs();
00173     const ID &theRetainedNodesID = theRetainedNodesDOFs->getID();    
00174 
00175     int size2 = RetainedDOFs.Size();
00176     for (int j=0; j<size2; j++) {
00177  int retained = RetainedDOFs(j);
00178  if (retained < 0 || retained >= theRetainedNode->getNumberDOF()) {
00179      cerr << "WARNING PenaltyMP_FE::setID(void) - unknown DOF ";
00180      cerr << retained << " at Node\n";
00181      myID(j+size1) = -1; // modify so nothing will be added
00182      result = -3;
00183  }     
00184  else {
00185      if (retained >= theRetainedNodesID.Size()) {
00186   cerr << "WARNING PenaltyMP_FE::setID(void) - ";
00187   cerr << " Nodes DOF_Group too small\n";
00188   myID(j+size1) = -1; // modify so nothing will be added 
00189   result = -4;
00190      }
00191      else
00192   myID(j+size1) = theRetainedNodesID(retained);
00193  }
00194     }
00195 
00196     myDOF_Groups(0) = theConstrainedNodesDOFs->getTag();
00197     myDOF_Groups(1) = theRetainedNodesDOFs->getTag();
00198 
00199     return result;
00200 }
00201 
00202 const Matrix &
00203 PenaltyMP_FE::getTangent(Integrator *theNewIntegrator)
00204 {
00205     if (theMP->isTimeVarying() == true)
00206  this->determineTangent();    
00207     return *tang;
00208 }
00209 
00210 const Vector &
00211 PenaltyMP_FE::getResidual(Integrator *theNewIntegrator)
00212 {
00213     // zero residual, CD = 0
00214     return *resid;
00215 }
00216 
00217 
00218 
00219 const Vector &
00220 PenaltyMP_FE::getTangForce(const Vector &disp, double fact)
00221 {
00222     // does nothing , zero residual for CD = 0
00223     return *resid;
00224 }
00225 
00226 void  
00227 PenaltyMP_FE::determineTangent(void)
00228 {
00229     // first determine [C] = [-I [Ccr]]
00230     C->Zero();
00231     const Matrix &constraint = theMP->getConstraint();
00232     int noRows = constraint.noRows();
00233     int noCols = constraint.noCols();
00234     
00235     for (int j=0; j<noRows; j++)
00236  (*C)(j,j) = -1.0;
00237     
00238     for (int i=0; i<noRows; i++)
00239  for (int j=0; j<noCols; j++)
00240      (*C)(i,j+noRows) = constraint(i,j);
00241     
00242 
00243     // now form the tangent: [K] = alpha * [C]^t[C]
00244     *(tang) = (*C)^(*C);
00245     *(tang) *= alpha;
00246 }
00247 
00248