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 #include "LOCA_Epetra_AnasaziOperator_Floquet.H"
00034 #include "Teuchos_ParameterList.hpp"
00035 #include "Teuchos_RCPDecl.hpp"
00036 #include "LOCA_GlobalData.H"
00037 #include "LOCA_ErrorCheck.H"
00038
00039 LOCA::Epetra::AnasaziOperator::Floquet::Floquet(
00040 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00041 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00042 const Teuchos::RCP<Teuchos::ParameterList>& eigenParams_,
00043 const Teuchos::RCP<Teuchos::ParameterList>& solverParams_,
00044 const Teuchos::RCP<NOX::Abstract::Group>& grp_)
00045 : globalData(global_data),
00046 myLabel("Floquet Transformation"),
00047 eigenParams(eigenParams_),
00048 solverParams(solverParams_),
00049 grp(grp_),
00050 xyztInterface()
00051 {
00052 string callingFunction =
00053 "LOCA::Epetra::AnasaziOperator::Floquet::Floquet()";
00054
00055 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00056 NOX::Abstract::Group::ReturnType status;
00057
00058 Teuchos::RCP<NOX::Epetra::Group> NEGrp
00059 = Teuchos::rcp_dynamic_cast<NOX::Epetra::Group>(grp);
00060 if (NEGrp == Teuchos::null) cout << callingFunction << " NEGrp cast failed " << endl;
00061 else cout << callingFunction << " NEGrp cast succeeded." << endl;
00062
00063 xyztInterface = Teuchos::rcp_dynamic_cast<LOCA::Epetra::Interface::xyzt>(NEGrp->getRequiredInterface());
00064 if (xyztInterface == Teuchos::null) cout << callingFunction << " xyztInterface cast failed " << endl;
00065 else cout << callingFunction << " xyztInterface cast succeeded." << endl;
00066
00067
00068
00069 xyztInterface->setFloquetFillFlag(true);
00070 status = grp->computeJacobian();
00071
00072 finalStatus =
00073 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00074 callingFunction);
00075 }
00076
00077 LOCA::Epetra::AnasaziOperator::Floquet::~Floquet()
00078 {
00079 xyztInterface->setFloquetFillFlag(false);
00080 }
00081
00082 const string&
00083 LOCA::Epetra::AnasaziOperator::Floquet::label() const
00084 {
00085 return myLabel;
00086 }
00087
00088 void
00089 LOCA::Epetra::AnasaziOperator::Floquet::apply(const NOX::Abstract::MultiVector& input,
00090 NOX::Abstract::MultiVector& output) const
00091 {
00092
00093
00094 Teuchos::RCP<NOX::Abstract::MultiVector> tmpVec = input.clone();
00095 for (int i=0; i < input.numVectors(); i++) {
00096 NOX::Abstract::Vector& nAV = tmpVec->operator[](i);
00097 NOX::Epetra::Vector& nEV = dynamic_cast<NOX::Epetra::Vector&>(nAV);
00098 Epetra_Vector& eV = nEV.getEpetraVector();
00099
00100 xyztInterface->beginFloquetOperatorApplication(eV);
00101 }
00102
00103
00104 NOX::Abstract::Group::ReturnType status =
00105 grp->applyJacobianInverseMultiVector(*solverParams, *(tmpVec.get()), output);
00106 globalData->locaErrorCheck->checkReturnType(status,
00107 "LOCA::Epetra::AnasaziOperator::Floquet::apply()");
00108
00109
00110 for (int i=0; i < input.numVectors(); i++) {
00111 NOX::Abstract::Vector& nAV = output.operator[](i);
00112 NOX::Epetra::Vector& nEV = dynamic_cast<NOX::Epetra::Vector&>(nAV);
00113 Epetra_Vector& eV = nEV.getEpetraVector();
00114
00115 xyztInterface->finishFloquetOperatorApplication(eV);
00116 }
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143 }
00144
00145 void
00146 LOCA::Epetra::AnasaziOperator::Floquet::transformEigenvalue(double& ev_r,
00147 double& ev_i) const
00148 {
00149
00150 }
00151
00152 NOX::Abstract::Group::ReturnType
00153 LOCA::Epetra::AnasaziOperator::Floquet::rayleighQuotient(
00154 const NOX::Abstract::Vector& evec_r,
00155 const NOX::Abstract::Vector& evec_i,
00156 double& rq_r, double& rq_i) const
00157 {
00158 string callingFunction =
00159 "LOCA::Epetra::AnasaziOperator::Floquet::rayleighQuotient()";
00160
00161
00162 Teuchos::RCP<NOX::Abstract::MultiVector> z = evec_r.createMultiVector(2, NOX::DeepCopy);
00163 z->operator[](1) = evec_i;
00164 Teuchos::RCP<NOX::Abstract::MultiVector> Az = z->clone(NOX::ShapeCopy);
00165
00166 apply(*(z.get()), *(Az.get()));
00167 NOX::Abstract::Vector& Ax = Az->operator[](0);
00168 NOX::Abstract::Vector& Ay = Az->operator[](1);
00169
00170 double mag = evec_r.innerProduct(evec_r) + evec_i.innerProduct(evec_i) ;
00171 rq_r = (Ax.innerProduct(evec_r) + Ay.innerProduct(evec_i)) / mag ;
00172 rq_i = (Ay.innerProduct(evec_r) - Ax.innerProduct(evec_i)) / mag ;
00173
00174 return NOX::Abstract::Group::Ok;
00175 }