EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraModelEval2DSim.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // EpetraExt: Epetra Extended - Linear Algebra Services Package
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #include "EpetraModelEval2DSim.hpp"
45 #include "Teuchos_ScalarTraits.hpp"
46 #include "Epetra_SerialComm.h"
47 #include "Epetra_CrsMatrix.h"
48 
50  const double d
51  ,const double p0
52  ,const double p1
53  ,const double x00
54  ,const double x01
55  ,const bool showGetInvalidArg
56  )
57  :d_(d),showGetInvalidArg_(showGetInvalidArg)
58 {
59  using Teuchos::rcp;
60 
61  epetra_comm_ = rcp(new Epetra_SerialComm());
62 
63  const int nx = 2;
64 
65  map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
66  x0_ = rcp(new Epetra_Vector(*map_x_)); (*x0_)[0] = x00; (*x0_)[1] = x01;
67  p_ = rcp(new Epetra_Vector(*map_x_)); (*p_)[0] = p0; (*p_)[1] = p1;
68 
69  // Initialize the graph for W CrsMatrix object
70  W_graph_ = rcp(new Epetra_CrsGraph(::Copy,*map_x_,nx));
71  {
72  int indices[nx] = { 0, 1 };
73  for( int i = 0; i < nx; ++i )
74  W_graph_->InsertGlobalIndices(i,nx,indices);
75  }
76  W_graph_->FillComplete();
77 
78  isInitialized_ = true;
79 
80 }
81 
82 // Overridden from EpetraExt::ModelEvaluator
83 
84 Teuchos::RCP<const Epetra_Map>
86 {
87  return map_x_;
88 }
89 
90 Teuchos::RCP<const Epetra_Map>
92 {
93  return map_x_;
94 }
95 
96 Teuchos::RCP<const Epetra_Vector>
98 {
99  return x0_;
100 }
101 
102 Teuchos::RCP<Epetra_Operator>
104 {
105  return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
106 }
107 
110 {
111  InArgsSetup inArgs;
112  inArgs.setModelEvalDescription(this->description());
113  inArgs.setSupports(IN_ARG_x,true);
114  return inArgs;
115 }
116 
119 {
120  OutArgsSetup outArgs;
121  outArgs.setModelEvalDescription(this->description());
122  outArgs.setSupports(OUT_ARG_f,true);
123  outArgs.setSupports(OUT_ARG_W,true);
124  outArgs.set_W_properties(
128  ,true // supportsAdjoint
129  )
130  );
131  return outArgs;
132 }
133 
134 void EpetraModelEval2DSim::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
135 {
136  using Teuchos::dyn_cast;
137  using Teuchos::rcp_dynamic_cast;
138  //
139  // Get the input arguments
140  //
141  const Epetra_Vector &x = *inArgs.get_x();
142  //
143  // Get the output arguments
144  //
145  Epetra_Vector *f_out = outArgs.get_f().get();
146  Epetra_Operator *W_out = outArgs.get_W().get();
147  if(showGetInvalidArg_) {
148  Epetra_Vector *g_out = outArgs.get_g(0).get();
149  }
150  //
151  // Compute the functions
152  //
153  const Epetra_Vector &p = *p_;
154  if(f_out) {
155  Epetra_Vector &f = *f_out;
156  f[0] = x[0] + x[1]*x[1] - p[0];
157  f[1] = d_ * ( x[0]*x[0] - x[1] - p[1] );
158  }
159  if(W_out) {
160  Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
161  DfDx.PutScalar(0.0);
162  //
163  // Fill W = DfDx
164  //
165  // W = DfDx = [ 1.0, 2*x[1] ]
166  // [ 2*d*x[0], -d ]
167  //
168  double values[2];
169  int indexes[2];
170  // Row [0]
171  values[0] = 1.0; indexes[0] = 0;
172  values[1] = 2.0*x[1]; indexes[1] = 1;
173  DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
174  // Row [1]
175  values[0] = 2.0*d_*x[0]; indexes[0] = 0;
176  values[1] = -d_; indexes[1] = 1;
177  DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
178  }
179 }
void set_W_properties(const DerivativeProperties &properties)
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 &lt;= j &amp;&amp; j &lt; this-&gt;Ng().
void setSupports(EOutArgsMembers arg, bool supports=true)
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
OutArgs createOutArgs() const
Teuchos::RCP< Epetra_Operator > get_W() const
void setSupports(EInArgsMembers arg, bool supports=true)
Teuchos::RCP< const Epetra_Map > get_f_map() const
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< const Epetra_Vector > get_x_init() const
EpetraModelEval2DSim(const double d=10.0, const double p0=2.0, const double p1=0.0, const double x00=1.0, const double x01=1.0, const bool showGetInvalidArg=false)
int PutScalar(double ScalarConstant)
void setModelEvalDescription(const std::string &modelEvalDescription)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Epetra_Map > get_x_map() const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluation< Epetra_Vector > get_f() const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.