00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_MultiContinuation_ArcLengthConstraint.H"
00043 #include "LOCA_MultiContinuation_ArcLengthGroup.H"
00044 #include "LOCA_GlobalData.H"
00045 #include "LOCA_ErrorCheck.H"
00046
00047 LOCA::MultiContinuation::ArcLengthConstraint::ArcLengthConstraint(
00048 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00049 const Teuchos::RCP<LOCA::MultiContinuation::ArcLengthGroup>& grp) :
00050 globalData(global_data),
00051 arcLengthGroup(grp),
00052 constraints(grp->getNumParams(), 1),
00053 isValidConstraints(false),
00054 conParamIDs(grp->getContinuationParameterIDs())
00055 {
00056 }
00057
00058 LOCA::MultiContinuation::ArcLengthConstraint::ArcLengthConstraint(
00059 const LOCA::MultiContinuation::ArcLengthConstraint& source,
00060 NOX::CopyType type) :
00061 globalData(source.globalData),
00062 arcLengthGroup(),
00063 constraints(source.constraints),
00064 isValidConstraints(false),
00065 conParamIDs(source.conParamIDs)
00066 {
00067 if (source.isValidConstraints && type == NOX::DeepCopy)
00068 isValidConstraints = true;
00069 }
00070
00071 LOCA::MultiContinuation::ArcLengthConstraint::~ArcLengthConstraint()
00072 {
00073 }
00074
00075 void
00076 LOCA::MultiContinuation::ArcLengthConstraint::setArcLengthGroup(const Teuchos::RCP<LOCA::MultiContinuation::ArcLengthGroup>& grp)
00077 {
00078 arcLengthGroup = grp;
00079 }
00080
00081 void
00082 LOCA::MultiContinuation::ArcLengthConstraint::copy(
00083 const LOCA::MultiContinuation::ConstraintInterface& src)
00084 {
00085 const LOCA::MultiContinuation::ArcLengthConstraint& source =
00086 dynamic_cast<const LOCA::MultiContinuation::ArcLengthConstraint&>(src);
00087
00088 if (this != &source) {
00089 globalData = source.globalData;
00090 constraints.assign(source.constraints);
00091 isValidConstraints = source.isValidConstraints;
00092 conParamIDs = source.conParamIDs;
00093 }
00094 }
00095
00096 Teuchos::RCP<LOCA::MultiContinuation::ConstraintInterface>
00097 LOCA::MultiContinuation::ArcLengthConstraint::clone(NOX::CopyType type) const
00098 {
00099 return Teuchos::rcp(new ArcLengthConstraint(*this, type));
00100 }
00101
00102 int
00103 LOCA::MultiContinuation::ArcLengthConstraint::numConstraints() const
00104 {
00105 return constraints.numRows();
00106 }
00107
00108 void
00109 LOCA::MultiContinuation::ArcLengthConstraint::setX(
00110 const NOX::Abstract::Vector& y)
00111 {
00112 isValidConstraints = false;
00113 }
00114
00115 void
00116 LOCA::MultiContinuation::ArcLengthConstraint::setParam(int paramID, double val)
00117 {
00118 isValidConstraints = false;
00119 }
00120
00121 void
00122 LOCA::MultiContinuation::ArcLengthConstraint::setParams(
00123 const vector<int>& paramIDs,
00124 const NOX::Abstract::MultiVector::DenseMatrix& vals)
00125 {
00126 isValidConstraints = false;
00127 }
00128
00129 NOX::Abstract::Group::ReturnType
00130 LOCA::MultiContinuation::ArcLengthConstraint::computeConstraints()
00131 {
00132 if (isValidConstraints)
00133 return NOX::Abstract::Group::Ok;
00134
00135 string callingFunction =
00136 "LOCA::MultiContinuation::ArcLengthConstraint::computeConstraints()";
00137 NOX::Abstract::Group::ReturnType status;
00138 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00139
00140
00141 if (!arcLengthGroup->isPredictor()) {
00142 status = arcLengthGroup->computePredictor();
00143 finalStatus =
00144 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00145 finalStatus,
00146 callingFunction);
00147 }
00148
00149
00150 const LOCA::MultiContinuation::ExtendedMultiVector& scaledTangent =
00151 arcLengthGroup->getScaledPredictorTangent();
00152 const LOCA::MultiContinuation::ExtendedMultiVector& tangent =
00153 arcLengthGroup->getPredictorTangent();
00154
00155
00156 Teuchos::RCP<LOCA::MultiContinuation::ExtendedMultiVector> secant =
00157 Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector>(tangent.clone(1));
00158 (*secant)[0].update(1.0, arcLengthGroup->getX(),
00159 -1.0, arcLengthGroup->getPrevX(), 0.0);
00160
00161
00162 secant->multiply(1.0, scaledTangent, constraints);
00163 for (int i=0; i<arcLengthGroup->getNumParams(); i++)
00164 constraints(i,0) -= arcLengthGroup->getStepSize(i) *
00165 scaledTangent[i].innerProduct(tangent[i]);
00166
00167 isValidConstraints = true;
00168
00169 return finalStatus;
00170 }
00171
00172 NOX::Abstract::Group::ReturnType
00173 LOCA::MultiContinuation::ArcLengthConstraint::computeDX()
00174 {
00175 if (!isValidConstraints)
00176 return computeConstraints();
00177 else
00178 return NOX::Abstract::Group::Ok;
00179 }
00180
00181 NOX::Abstract::Group::ReturnType
00182 LOCA::MultiContinuation::ArcLengthConstraint::computeDP(
00183 const vector<int>& paramIDs,
00184 NOX::Abstract::MultiVector::DenseMatrix& dgdp,
00185 bool isValidG)
00186 {
00187 string callingFunction =
00188 "LOCA::MultiContinuation::ArcLengthConstraint::computeDP()";
00189 NOX::Abstract::Group::ReturnType status;
00190 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00191
00192
00193 if (!isValidG && !isValidConstraints) {
00194 status = computeConstraints();
00195 finalStatus =
00196 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00197 finalStatus,
00198 callingFunction);
00199 }
00200 if (!isValidG) {
00201 for (int i=0; i<constraints.numRows(); i++)
00202 dgdp(i,0) = constraints(i,0);
00203 }
00204
00205
00206 const LOCA::MultiContinuation::ExtendedMultiVector& scaledTangent =
00207 arcLengthGroup->getScaledPredictorTangent();
00208
00209
00210
00211
00212 vector<int>::const_iterator it;
00213 int idx;
00214 for (unsigned int i=0; i<paramIDs.size(); i++) {
00215 it = find(conParamIDs.begin(), conParamIDs.end(), paramIDs[i]);
00216 if (it == conParamIDs.end())
00217 for (int k=0; k<constraints.numRows(); k++)
00218 dgdp(k,i+1) = 0.0;
00219 else {
00220 idx = it - conParamIDs.begin();
00221 for (int k=0; k<constraints.numRows(); k++)
00222 dgdp(k,i+1) = scaledTangent.getScalar(k,idx);
00223 }
00224 }
00225
00226 return finalStatus;
00227 }
00228
00229 bool
00230 LOCA::MultiContinuation::ArcLengthConstraint::isConstraints() const
00231 {
00232 return isValidConstraints;
00233 }
00234
00235 bool
00236 LOCA::MultiContinuation::ArcLengthConstraint::isDX() const
00237 {
00238 return isValidConstraints;
00239 }
00240
00241 const NOX::Abstract::MultiVector::DenseMatrix&
00242 LOCA::MultiContinuation::ArcLengthConstraint::getConstraints() const
00243 {
00244 return constraints;
00245 }
00246
00247 const NOX::Abstract::MultiVector*
00248 LOCA::MultiContinuation::ArcLengthConstraint::getDX() const
00249 {
00250
00251 const LOCA::MultiContinuation::ExtendedMultiVector& tangent =
00252 arcLengthGroup->getScaledPredictorTangent();
00253
00254 return tangent.getXMultiVec().get();
00255 }
00256
00257 bool
00258 LOCA::MultiContinuation::ArcLengthConstraint::isDXZero() const
00259 {
00260 return false;
00261 }