EpetraModelEval2DSim.cpp

00001 #include "EpetraModelEval2DSim.hpp"
00002 #include "Teuchos_ScalarTraits.hpp"
00003 #include "Epetra_SerialComm.h"
00004 #include "Epetra_CrsMatrix.h"
00005 
00006 EpetraModelEval2DSim::EpetraModelEval2DSim(
00007   const double         d
00008   ,const double        p0
00009   ,const double        p1
00010   ,const double        x00
00011   ,const double        x01
00012   ,const bool          showGetInvalidArg
00013   )
00014   :d_(d),showGetInvalidArg_(showGetInvalidArg)
00015 {
00016   using Teuchos::rcp;
00017 
00018   epetra_comm_ = rcp(new Epetra_SerialComm());
00019 
00020   const int nx = 2;
00021 
00022   map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
00023   x0_ = rcp(new Epetra_Vector(*map_x_));  (*x0_)[0] = x00; (*x0_)[1] = x01;
00024   p_  = rcp(new Epetra_Vector(*map_x_));  (*p_)[0] = p0; (*p_)[1] = p1;
00025 
00026   // Initialize the graph for W CrsMatrix object
00027   W_graph_ = rcp(new Epetra_CrsGraph(::Copy,*map_x_,nx));
00028   {
00029     int indices[nx] = { 0, 1 };
00030     for( int i = 0; i < nx; ++i )
00031       W_graph_->InsertGlobalIndices(i,nx,indices);
00032   }
00033   W_graph_->FillComplete();
00034 
00035   isInitialized_ = true;
00036 
00037 }
00038 
00039 // Overridden from EpetraExt::ModelEvaluator
00040 
00041 Teuchos::RefCountPtr<const Epetra_Map>
00042 EpetraModelEval2DSim::get_x_map() const
00043 {
00044   return map_x_;
00045 }
00046 
00047 Teuchos::RefCountPtr<const Epetra_Map>
00048 EpetraModelEval2DSim::get_f_map() const
00049 {
00050   return map_x_;
00051 }
00052 
00053 Teuchos::RefCountPtr<const Epetra_Vector>
00054 EpetraModelEval2DSim::get_x_init() const
00055 {
00056   return x0_;
00057 }
00058 
00059 Teuchos::RefCountPtr<Epetra_Operator>
00060 EpetraModelEval2DSim::create_W() const
00061 {
00062   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00063 }
00064 
00065 EpetraExt::ModelEvaluator::InArgs
00066 EpetraModelEval2DSim::createInArgs() const
00067 {
00068   InArgsSetup inArgs;
00069   inArgs.setModelEvalDescription(this->description());
00070   inArgs.setSupports(IN_ARG_x,true);
00071   return inArgs;
00072 }
00073 
00074 EpetraExt::ModelEvaluator::OutArgs
00075 EpetraModelEval2DSim::createOutArgs() const
00076 {
00077   OutArgsSetup outArgs;
00078   outArgs.setModelEvalDescription(this->description());
00079   outArgs.setSupports(OUT_ARG_f,true);
00080   outArgs.setSupports(OUT_ARG_W,true);
00081   outArgs.set_W_properties(
00082     DerivativeProperties(
00083       DERIV_LINEARITY_NONCONST
00084       ,DERIV_RANK_FULL
00085       ,true // supportsAdjoint
00086       )
00087     );
00088   return outArgs;
00089 }
00090 
00091 void EpetraModelEval2DSim::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
00092 {
00093   using Teuchos::dyn_cast;
00094   using Teuchos::rcp_dynamic_cast;
00095   //
00096   // Get the input arguments
00097   //
00098   const Epetra_Vector &x = *inArgs.get_x();
00099   //
00100   // Get the output arguments
00101   //
00102   Epetra_Vector       *f_out = outArgs.get_f().get();
00103   Epetra_Operator     *W_out = outArgs.get_W().get();
00104   if(showGetInvalidArg_) {
00105     Epetra_Vector *g_out = outArgs.get_g(0).get();
00106   }
00107   //
00108   // Compute the functions
00109   //
00110   const Epetra_Vector &p = *p_;
00111   if(f_out) {
00112     Epetra_Vector &f = *f_out;
00113     f[0] =        x[0]      + x[1]*x[1] - p[0];
00114     f[1] = d_ * ( x[0]*x[0] - x[1]      - p[1] );
00115   }
00116   if(W_out) {
00117     Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
00118     DfDx.PutScalar(0.0);
00119     //
00120     // Fill W = DfDx
00121     //
00122     // W = DfDx = [      1.0,  2*x[1] ]
00123     //            [ 2*d*x[0],     -d  ]
00124     //
00125     double values[2];
00126     int indexes[2];
00127     // Row [0]
00128     values[0] = 1.0;           indexes[0] = 0;
00129     values[1] = 2.0*x[1];      indexes[1] = 1;
00130     DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
00131     // Row [1]
00132     values[0] = 2.0*d_*x[0];   indexes[0] = 0;
00133     values[1] = -d_;           indexes[1] = 1;
00134     DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
00135   }
00136 }

Generated on Thu Feb 11 16:28:43 2010 for EpetraExt/Thyra Adapters by  doxygen 1.4.7