Thyra Package Browser (Single Doxygen Collection) Version of the Day
Simple2DTpetraModelEvaluator_def.hpp
Go to the documentation of this file.
00001 
00002 #ifndef SIMPLE_2D_TPETRA_MODEL_EVALUATOR_HPP
00003 #define SIMPLE_2D_TPETRA_MODEL_EVALUATOR_HPP
00004 
00005 
00006 #include "Simple2DTpetraModelEvaluator_decl.hpp"
00007 #include "Thyra_TpetraThyraWrappers.hpp"
00008 #include "Tpetra_CrsMatrix.hpp"
00009 #include "Teuchos_DefaultComm.hpp"
00010 
00011 
00012 // Constructors/Initializers/Accessors
00013 
00014 
00015 template<class Scalar>
00016 Simple2DTpetraModelEvaluator<Scalar>::Simple2DTpetraModelEvaluator()
00017   : d_(0.0)
00018 {
00019 
00020   using Teuchos::RCP;
00021   using Teuchos::rcp;
00022   using Teuchos::tuple;
00023   using Thyra::VectorBase;
00024   using Thyra::createMember;
00025   typedef Thyra::ModelEvaluatorBase MEB;
00026   typedef Teuchos::ScalarTraits<Scalar> ST;
00027 
00028   const int dim = 2;
00029 
00030   //
00031   // A) Create the structure for the problem
00032   //
00033   
00034   MEB::InArgsSetup<Scalar> inArgs;
00035   inArgs.setModelEvalDescription(this->description());
00036   inArgs.setSupports(MEB::IN_ARG_x);
00037   prototypeInArgs_ = inArgs;
00038   
00039   MEB::OutArgsSetup<Scalar> outArgs;
00040   outArgs.setModelEvalDescription(this->description());
00041   outArgs.setSupports(MEB::OUT_ARG_f);
00042   outArgs.setSupports(MEB::OUT_ARG_W_op);
00043   prototypeOutArgs_ = outArgs;
00044 
00045   //
00046   // B) Create the Tpetra objects
00047   //
00048 
00049   const RCP<const Teuchos::Comm<int> > comm =
00050     Teuchos::DefaultComm<int>::getComm();
00051   
00052   const RCP<const Tpetra::Map<int> > map = 
00053     rcp(new Tpetra::Map<int>(dim, 0, comm));
00054 
00055   W_op_graph_ = rcp(new Tpetra::CrsGraph<int>(map, dim));
00056   W_op_graph_->insertGlobalIndices(0, tuple<int>(0, 1)());
00057   W_op_graph_->insertGlobalIndices(1, tuple<int>(0, 1)());
00058   W_op_graph_->fillComplete();
00059 
00060   p_.resize(dim, ST::zero());
00061 
00062   x0_ = rcp(new Tpetra::Vector<Scalar, int>(map));
00063   x0_->putScalar(ST::zero());
00064 
00065   //
00066   // C) Create the Thyra wrapped objects
00067   //
00068 
00069   x_space_ = Thyra::createVectorSpace<Scalar>(map);
00070   f_space_ = x_space_;
00071 
00072   nominalValues_ = inArgs;
00073   nominalValues_.set_x(Thyra::createVector(x0_, x_space_));
00074 
00075   //
00076   // D) Set initial values through interface functions
00077   //
00078 
00079   set_d(10.0);
00080   set_p(Teuchos::tuple<Scalar>(2.0, 0.0)());
00081   set_x0(Teuchos::tuple<Scalar>(1.0, 1.0)());
00082 
00083 }
00084 
00085 
00086 template<class Scalar>
00087 void Simple2DTpetraModelEvaluator<Scalar>::set_d(const Scalar &d)
00088 {
00089   d_ = d;
00090 }
00091 
00092 
00093 template<class Scalar>
00094 void Simple2DTpetraModelEvaluator<Scalar>::set_p(const Teuchos::ArrayView<const Scalar> &p)
00095 {
00096 #ifdef TEUCHOS_DEBUG
00097   TEUCHOS_ASSERT_EQUALITY(p_.size(), p.size());
00098 #endif
00099   p_().assign(p);
00100 }
00101 
00102 
00103 template<class Scalar>
00104 void Simple2DTpetraModelEvaluator<Scalar>::set_x0(const Teuchos::ArrayView<const Scalar> &x0_in)
00105 {
00106 #ifdef TEUCHOS_DEBUG
00107   TEUCHOS_ASSERT_EQUALITY(x_space_->dim(), x0_in.size());
00108 #endif
00109   x0_->get1dViewNonConst()().assign(x0_in);
00110 }
00111 
00112 
00113 // Public functions overridden from ModelEvaulator
00114 
00115 
00116 template<class Scalar>
00117 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00118 Simple2DTpetraModelEvaluator<Scalar>::get_x_space() const
00119 {
00120   return x_space_;
00121 }
00122 
00123 
00124 template<class Scalar>
00125 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00126 Simple2DTpetraModelEvaluator<Scalar>::get_f_space() const
00127 {
00128   return f_space_;
00129 }
00130 
00131 
00132 template<class Scalar>
00133 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00134 Simple2DTpetraModelEvaluator<Scalar>::getNominalValues() const
00135 {
00136   return nominalValues_;
00137 }
00138 
00139 
00140 template<class Scalar>
00141 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00142 Simple2DTpetraModelEvaluator<Scalar>::create_W_op() const
00143 {
00144   return Thyra::createLinearOp(
00145     Teuchos::RCP<Tpetra::Operator<Scalar,int> >(
00146       Teuchos::rcp(new Tpetra::CrsMatrix<Scalar,int>(W_op_graph_))
00147       )
00148     );
00149 }
00150 
00151 
00152 template<class Scalar>
00153 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00154 Simple2DTpetraModelEvaluator<Scalar>::createInArgs() const
00155 {
00156   return prototypeInArgs_;
00157 }
00158 
00159 
00160 // Private functions overridden from ModelEvaulatorDefaultBase
00161 
00162 
00163 template<class Scalar>
00164 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00165 Simple2DTpetraModelEvaluator<Scalar>::createOutArgsImpl() const
00166 {
00167   return prototypeOutArgs_;
00168 }
00169 
00170 
00171 template<class Scalar>
00172 void Simple2DTpetraModelEvaluator<Scalar>::evalModelImpl(
00173   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00174   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00175   ) const
00176 {
00177   using Teuchos::RCP;
00178   using Teuchos::ArrayRCP;
00179   using Teuchos::Array;
00180   using Teuchos::tuple;
00181   using Teuchos::rcp_dynamic_cast;
00182   typedef Teuchos::ScalarTraits<Scalar> ST;
00183   typedef Thyra::TpetraOperatorVectorExtraction<Scalar,int> ConverterT;
00184 
00185   const RCP<const Tpetra::Vector<Scalar, int> > x_vec =
00186     ConverterT::getConstTpetraVector(inArgs.get_x());
00187   const ArrayRCP<const Scalar> x = x_vec->get1dView();
00188 
00189   const RCP<Tpetra::Vector<Scalar, int> > f_vec =
00190     ConverterT::getTpetraVector(outArgs.get_f());
00191 
00192   const RCP<Tpetra::CrsMatrix<Scalar, int> > W =
00193     rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,int> >(
00194       ConverterT::getTpetraOperator(outArgs.get_W_op()),
00195       true
00196       );
00197 
00198   if (nonnull(f_vec)) {
00199     const ArrayRCP<Scalar> f = f_vec->get1dViewNonConst();
00200     f[0] = x[0] + x[1]*x[1] - p_[0];
00201     f[1] = d_ * (x[0]*x[0] -x[1] - p_[1]);
00202   }
00203 
00204   if (nonnull(W)) {
00205     W->setAllToScalar(ST::zero());
00206     W->sumIntoGlobalValues(0, tuple<int>(0, 1), tuple<Scalar>(1.0, 2.0*x[1]));
00207     W->sumIntoGlobalValues(1, tuple<int>(0, 1), tuple<Scalar>(2.0*d_*x[0], -d_));
00208   }
00209 
00210 }
00211 
00212 
00213 #endif // SIMPLE_2D_TPETRA_MODEL_EVALUATOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines