Thyra Package Browser (Single Doxygen Collection) Version of the Day
Simple2DTpetraModelEvaluator_def.hpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 
00045 #ifndef SIMPLE_2D_TPETRA_MODEL_EVALUATOR_HPP
00046 #define SIMPLE_2D_TPETRA_MODEL_EVALUATOR_HPP
00047 
00048 
00049 #include "Simple2DTpetraModelEvaluator_decl.hpp"
00050 #include "Thyra_TpetraThyraWrappers.hpp"
00051 #include "Tpetra_CrsMatrix.hpp"
00052 #include "Teuchos_DefaultComm.hpp"
00053 
00054 
00055 // Constructors/Initializers/Accessors
00056 
00057 
00058 template<class Scalar>
00059 Simple2DTpetraModelEvaluator<Scalar>::Simple2DTpetraModelEvaluator()
00060   : d_(0.0)
00061 {
00062 
00063   using Teuchos::RCP;
00064   using Teuchos::rcp;
00065   using Teuchos::tuple;
00066   using Thyra::VectorBase;
00067   using Thyra::createMember;
00068   typedef Thyra::ModelEvaluatorBase MEB;
00069   typedef Teuchos::ScalarTraits<Scalar> ST;
00070 
00071   const int dim = 2;
00072 
00073   //
00074   // A) Create the structure for the problem
00075   //
00076   
00077   MEB::InArgsSetup<Scalar> inArgs;
00078   inArgs.setModelEvalDescription(this->description());
00079   inArgs.setSupports(MEB::IN_ARG_x);
00080   prototypeInArgs_ = inArgs;
00081   
00082   MEB::OutArgsSetup<Scalar> outArgs;
00083   outArgs.setModelEvalDescription(this->description());
00084   outArgs.setSupports(MEB::OUT_ARG_f);
00085   outArgs.setSupports(MEB::OUT_ARG_W_op);
00086   prototypeOutArgs_ = outArgs;
00087 
00088   //
00089   // B) Create the Tpetra objects
00090   //
00091 
00092   const RCP<const Teuchos::Comm<int> > comm =
00093     Teuchos::DefaultComm<int>::getComm();
00094   
00095   const RCP<const Tpetra::Map<int> > map = 
00096     rcp(new Tpetra::Map<int>(dim, 0, comm));
00097 
00098   W_op_graph_ = rcp(new Tpetra::CrsGraph<int>(map, dim));
00099   W_op_graph_->insertGlobalIndices(0, tuple<int>(0, 1)());
00100   W_op_graph_->insertGlobalIndices(1, tuple<int>(0, 1)());
00101   W_op_graph_->fillComplete();
00102 
00103   p_.resize(dim, ST::zero());
00104 
00105   x0_ = rcp(new Tpetra::Vector<Scalar, int>(map));
00106   x0_->putScalar(ST::zero());
00107 
00108   //
00109   // C) Create the Thyra wrapped objects
00110   //
00111 
00112   x_space_ = Thyra::createVectorSpace<Scalar>(map);
00113   f_space_ = x_space_;
00114 
00115   nominalValues_ = inArgs;
00116   nominalValues_.set_x(Thyra::createVector(x0_, x_space_));
00117 
00118   //
00119   // D) Set initial values through interface functions
00120   //
00121 
00122   set_d(10.0);
00123   set_p(Teuchos::tuple<Scalar>(2.0, 0.0)());
00124   set_x0(Teuchos::tuple<Scalar>(1.0, 1.0)());
00125 
00126 }
00127 
00128 
00129 template<class Scalar>
00130 void Simple2DTpetraModelEvaluator<Scalar>::set_d(const Scalar &d)
00131 {
00132   d_ = d;
00133 }
00134 
00135 
00136 template<class Scalar>
00137 void Simple2DTpetraModelEvaluator<Scalar>::set_p(const Teuchos::ArrayView<const Scalar> &p)
00138 {
00139 #ifdef TEUCHOS_DEBUG
00140   TEUCHOS_ASSERT_EQUALITY(p_.size(), p.size());
00141 #endif
00142   p_().assign(p);
00143 }
00144 
00145 
00146 template<class Scalar>
00147 void Simple2DTpetraModelEvaluator<Scalar>::set_x0(const Teuchos::ArrayView<const Scalar> &x0_in)
00148 {
00149 #ifdef TEUCHOS_DEBUG
00150   TEUCHOS_ASSERT_EQUALITY(x_space_->dim(), x0_in.size());
00151 #endif
00152   x0_->get1dViewNonConst()().assign(x0_in);
00153 }
00154 
00155 
00156 // Public functions overridden from ModelEvaulator
00157 
00158 
00159 template<class Scalar>
00160 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00161 Simple2DTpetraModelEvaluator<Scalar>::get_x_space() const
00162 {
00163   return x_space_;
00164 }
00165 
00166 
00167 template<class Scalar>
00168 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00169 Simple2DTpetraModelEvaluator<Scalar>::get_f_space() const
00170 {
00171   return f_space_;
00172 }
00173 
00174 
00175 template<class Scalar>
00176 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00177 Simple2DTpetraModelEvaluator<Scalar>::getNominalValues() const
00178 {
00179   return nominalValues_;
00180 }
00181 
00182 
00183 template<class Scalar>
00184 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00185 Simple2DTpetraModelEvaluator<Scalar>::create_W_op() const
00186 {
00187   return Thyra::createLinearOp(
00188     Teuchos::RCP<Tpetra::Operator<Scalar,int> >(
00189       Teuchos::rcp(new Tpetra::CrsMatrix<Scalar,int>(W_op_graph_))
00190       )
00191     );
00192 }
00193 
00194 
00195 template<class Scalar>
00196 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00197 Simple2DTpetraModelEvaluator<Scalar>::createInArgs() const
00198 {
00199   return prototypeInArgs_;
00200 }
00201 
00202 
00203 // Private functions overridden from ModelEvaulatorDefaultBase
00204 
00205 
00206 template<class Scalar>
00207 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00208 Simple2DTpetraModelEvaluator<Scalar>::createOutArgsImpl() const
00209 {
00210   return prototypeOutArgs_;
00211 }
00212 
00213 
00214 template<class Scalar>
00215 void Simple2DTpetraModelEvaluator<Scalar>::evalModelImpl(
00216   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00217   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00218   ) const
00219 {
00220   using Teuchos::RCP;
00221   using Teuchos::ArrayRCP;
00222   using Teuchos::Array;
00223   using Teuchos::tuple;
00224   using Teuchos::rcp_dynamic_cast;
00225   typedef Teuchos::ScalarTraits<Scalar> ST;
00226   typedef Thyra::TpetraOperatorVectorExtraction<Scalar,int> ConverterT;
00227 
00228   const RCP<const Tpetra::Vector<Scalar, int> > x_vec =
00229     ConverterT::getConstTpetraVector(inArgs.get_x());
00230   const ArrayRCP<const Scalar> x = x_vec->get1dView();
00231 
00232   const RCP<Tpetra::Vector<Scalar, int> > f_vec =
00233     ConverterT::getTpetraVector(outArgs.get_f());
00234 
00235   const RCP<Tpetra::CrsMatrix<Scalar, int> > W =
00236     rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,int> >(
00237       ConverterT::getTpetraOperator(outArgs.get_W_op()),
00238       true
00239       );
00240 
00241   if (nonnull(f_vec)) {
00242     const ArrayRCP<Scalar> f = f_vec->get1dViewNonConst();
00243     f[0] = x[0] + x[1]*x[1] - p_[0];
00244     f[1] = d_ * (x[0]*x[0] -x[1] - p_[1]);
00245   }
00246 
00247   if (nonnull(W)) {
00248     W->setAllToScalar(ST::zero());
00249     W->sumIntoGlobalValues(0, tuple<int>(0, 1), tuple<Scalar>(1.0, 2.0*x[1]));
00250     W->sumIntoGlobalValues(1, tuple<int>(0, 1), tuple<Scalar>(2.0*d_*x[0], -d_));
00251   }
00252 
00253 }
00254 
00255 
00256 #endif // SIMPLE_2D_TPETRA_MODEL_EVALUATOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines