Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_SGQuadModelEvaluator.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
44 #include "Stokhos_Quadrature.hpp"
48 #include "Epetra_Map.h"
49 #include "Epetra_Vector.h"
50 #include "Teuchos_TimeMonitor.hpp"
51 #include "Teuchos_Assert.hpp"
52 
56  me(me_),
57  num_p(0),
58  num_g(0),
59  x_dot_qp(),
60  x_qp(),
61  p_qp(),
62  f_qp(),
63  W_qp(),
64  dfdp_qp(),
65  g_qp(),
66  dgdx_qp(),
67  dgdx_dot_qp(),
68  dgdp_qp()
69 {
70  // Create storage for x_dot, x, and p at a quad point
71  InArgs me_inargs = me->createInArgs();
72  num_p = me_inargs.Np();
73  if (me_inargs.supports(IN_ARG_x_dot))
74  x_dot_qp = Teuchos::rcp(new Epetra_Vector(*(me->get_x_map())));
75  if (me_inargs.supports(IN_ARG_x))
76  x_qp = Teuchos::rcp(new Epetra_Vector((*me->get_x_map())));
77  p_qp.resize(num_p);
78  for (int i=0; i<num_p; i++)
79  p_qp[i] = Teuchos::rcp(new Epetra_Vector(*(me->get_p_map(i))));
80 
81  // Create storage for f and W at a quad point
82  OutArgs me_outargs = me->createOutArgs();
83  num_g = me_outargs.Ng();
84 
85  // f
86  if (me_outargs.supports(OUT_ARG_f))
87  f_qp = Teuchos::rcp(new Epetra_Vector(*(me->get_f_map())));
88 
89  // W
90  if (me_outargs.supports(OUT_ARG_W))
91  W_qp = me->create_W();
92 
93  // df/dp
94  dfdp_qp.resize(num_p);
95  for (int i=0; i<num_p; i++)
96  if (me_outargs.supports(OUT_ARG_DfDp,i).supports(DERIV_MV_BY_COL))
97  dfdp_qp[i] = EpetraExt::ModelEvaluator::Derivative(
99  *(me->get_f_map()),
100  me->get_p_map(i)->NumGlobalElements())));
101  else if (me_outargs.supports(OUT_ARG_DfDp,i).supports(DERIV_TRANS_MV_BY_ROW))
102  dfdp_qp[i] = EpetraExt::ModelEvaluator::Derivative(
104  *(me->get_p_map(i)),
105  me->get_f_map()->NumGlobalElements())));
106  else if (me_outargs.supports(OUT_ARG_DfDp,i).supports(DERIV_LINEAR_OP))
107  dfdp_qp[i] = EpetraExt::ModelEvaluator::Derivative(
108  me->create_DfDp_op(i));
109 
110  g_qp.resize(num_g);
113  dgdp_qp.resize(num_g);
114  for (int i=0; i<num_g; i++) {
115 
116  // g
117  g_qp[i] =
118  Teuchos::rcp(new Epetra_Vector(*(me->get_g_map(i))));
119 
120  // dg/dx
121  if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_TRANS_MV_BY_ROW))
122  dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
124  *(me->get_x_map()),
125  me->get_g_map(i)->NumGlobalElements())));
126  else if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_MV_BY_COL))
127  dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
129  *(me->get_g_map(i)),
130  me->get_x_map()->NumGlobalElements())));
131  else if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP))
132  dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
133  me->create_DgDx_op(i));
134 
135  // dg/dx_dot
136  if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_TRANS_MV_BY_ROW))
137  dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
139  *(me->get_x_map()),
140  me->get_g_map(i)->NumGlobalElements())));
141  else if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_MV_BY_COL))
142  dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
144  *(me->get_g_map(i)),
145  me->get_x_map()->NumGlobalElements())));
146  else if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP))
147  dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
148  me->create_DgDx_dot_op(i));
149 
150  // dg/dp
151  dgdp_qp[i].resize(num_p);
152  for (int j=0; j<num_p; j++)
153  if (me_outargs.supports(OUT_ARG_DgDp, i, j).supports(DERIV_TRANS_MV_BY_ROW))
154  dgdp_qp[i][j] = EpetraExt::ModelEvaluator::Derivative(
156  *(me->get_p_map(j)),
157  me->get_g_map(i)->NumGlobalElements())));
158  else if (me_outargs.supports(OUT_ARG_DgDp, i, j).supports(DERIV_MV_BY_COL))
159  dgdp_qp[i][j] = EpetraExt::ModelEvaluator::Derivative(
161  *(me->get_g_map(i)),
162  me->get_p_map(j)->NumGlobalElements())));
163  else if (me_outargs.supports(OUT_ARG_DgDp, i, j).supports(DERIV_LINEAR_OP))
164  dgdp_qp[i][j] = EpetraExt::ModelEvaluator::Derivative(
165  me->create_DgDp_op(i,j));
166  }
167 }
168 
169 // Overridden from EpetraExt::ModelEvaluator
170 
173 get_x_map() const
174 {
175  return me->get_x_map();
176 }
177 
180 get_f_map() const
181 {
182  return me->get_f_map();
183 }
184 
187 get_p_map(int l) const
188 {
189  return me->get_p_map(l);
190 }
191 
194 get_g_map(int l) const
195 {
196  return me->get_g_map(l);
197 }
198 
201 get_p_names(int l) const
202 {
203  return me->get_p_names(l);
204 }
205 
208 get_x_init() const
209 {
210  return me->get_x_init();
211 }
212 
215 get_p_init(int l) const
216 {
217  return me->get_p_init(l);
218 }
219 
222 create_W() const
223 {
224  return me->create_W();
225 }
226 
227 EpetraExt::ModelEvaluator::InArgs
230 {
231  InArgsSetup inArgs;
232  InArgs me_inargs = me->createInArgs();
233 
234  inArgs.setModelEvalDescription(this->description());
235  inArgs.set_Np(num_p);
236  inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot));
237  inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x));
238  inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
239  inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
240  inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
241 
242  for (int i=0; i<num_p; i++)
243  inArgs.setSupports(IN_ARG_p_sg, i, true);
244  inArgs.setSupports(IN_ARG_x_sg, me_inargs.supports(IN_ARG_x));
245  inArgs.setSupports(IN_ARG_x_dot_sg, me_inargs.supports(IN_ARG_x_dot));
246  inArgs.setSupports(IN_ARG_sg_basis, true);
247  inArgs.setSupports(IN_ARG_sg_quadrature, true);
248 
249  return inArgs;
250 }
251 
252 EpetraExt::ModelEvaluator::OutArgs
255 {
256  OutArgsSetup outArgs;
257  OutArgs me_outargs = me->createOutArgs();
258 
259  outArgs.setModelEvalDescription(this->description());
260  outArgs.set_Np_Ng(num_p, num_g);
261  outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f));
262  outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W));
263  for (int j=0; j<num_p; j++)
264  outArgs.setSupports(OUT_ARG_DfDp, j,
265  me_outargs.supports(OUT_ARG_DfDp, j));
266  for (int i=0; i<num_g; i++) {
267  outArgs.setSupports(OUT_ARG_DgDx, i,
268  me_outargs.supports(OUT_ARG_DgDx, i));
269  outArgs.setSupports(OUT_ARG_DgDx_dot, i,
270  me_outargs.supports(OUT_ARG_DgDx_dot, i));
271  for (int j=0; j<num_p; j++)
272  outArgs.setSupports(OUT_ARG_DgDp, i, j,
273  me_outargs.supports(OUT_ARG_DgDp, i, j));
274  }
275 
276  outArgs.setSupports(OUT_ARG_f_sg, me_outargs.supports(OUT_ARG_f));
277  if (me_outargs.supports(OUT_ARG_W)) {
278  outArgs.set_W_properties(me_outargs.get_W_properties());
279  outArgs.setSupports(OUT_ARG_W_sg, true);
280  }
281  for (int j=0; j<num_p; j++) {
282  outArgs.setSupports(OUT_ARG_DfDp_sg, j,
283  me_outargs.supports(OUT_ARG_DfDp, j));
284  }
285  for (int i=0; i<num_g; i++) {
286  outArgs.setSupports(OUT_ARG_g_sg, i, true);
287  outArgs.setSupports(OUT_ARG_DgDx_sg, i,
288  me_outargs.supports(OUT_ARG_DgDx, i));
289  outArgs.setSupports(OUT_ARG_DgDx_dot_sg, i,
290  me_outargs.supports(OUT_ARG_DgDx_dot, i));
291  for (int j=0; j<num_p; j++)
292  outArgs.setSupports(OUT_ARG_DgDp_sg, i, j,
293  me_outargs.supports(OUT_ARG_DgDp, i, j));
294  }
295 
296  return outArgs;
297 }
298 
299 void
301 evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
302 {
303  // Create underlying inargs
304  InArgs me_inargs = me->createInArgs();
305  if (me_inargs.supports(IN_ARG_x))
306  me_inargs.set_x(inArgs.get_x());
307  if (me_inargs.supports(IN_ARG_x_dot))
308  me_inargs.set_x_dot(inArgs.get_x_dot());
309  if (me_inargs.supports(IN_ARG_alpha))
310  me_inargs.set_alpha(inArgs.get_alpha());
311  if (me_inargs.supports(IN_ARG_beta))
312  me_inargs.set_beta(inArgs.get_beta());
313  if (me_inargs.supports(IN_ARG_t))
314  me_inargs.set_t(inArgs.get_t());
315  for (int i=0; i<num_p; i++)
316  me_inargs.set_p(i, inArgs.get_p(i));
317 
318  // Create underlying outargs
319  OutArgs me_outargs = me->createOutArgs();
320  if (me_outargs.supports(OUT_ARG_f))
321  me_outargs.set_f(outArgs.get_f());
322  if (me_outargs.supports(OUT_ARG_W))
323  me_outargs.set_W(outArgs.get_W());
324  for (int j=0; j<num_p; j++)
325  if (!outArgs.supports(OUT_ARG_DfDp, j).none())
326  me_outargs.set_DfDp(j, outArgs.get_DfDp(j));
327  for (int i=0; i<num_g; i++) {
328  me_outargs.set_g(i, outArgs.get_g(i));
329  if (!outArgs.supports(OUT_ARG_DgDx, i).none())
330  me_outargs.set_DgDx(i, outArgs.get_DgDx(i));
331  if (!outArgs.supports(OUT_ARG_DgDx_dot, i).none())
332  me_outargs.set_DgDx(i, outArgs.get_DgDx_dot(i));
333  for (int j=0; j<num_p; j++)
334  if (!outArgs.supports(OUT_ARG_DgDp, i, j).none())
335  me_outargs.set_DgDp(i, j, outArgs.get_DgDp(i,j));
336  }
337 
338  bool do_quad = false;
339  InArgs::sg_const_vector_t x_sg;
340  InArgs::sg_const_vector_t x_dot_sg;
342  OutArgs::sg_vector_t f_sg;
343  OutArgs::sg_operator_t W_sg;
344  Teuchos::Array<SGDerivative> dfdp_sg(num_p);
346  Teuchos::Array<SGDerivative> dgdx_sg(num_g);
347  Teuchos::Array<SGDerivative> dgdx_dot_sg(num_g);
349  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_basis() == Teuchos::null,
350  std::logic_error,
351  "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
352  "SG basis inArg cannot be null!");
353  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_quadrature() == Teuchos::null,
354  std::logic_error,
355  "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
356  "SG quadrature inArg cannot be null!");
358  inArgs.get_sg_basis();
360  inArgs.get_sg_quadrature();
361  if (inArgs.supports(IN_ARG_x_sg)) {
362  x_sg = inArgs.get_x_sg();
363  if (x_sg != Teuchos::null) {
364  do_quad = true;
365  }
366  }
367  if (inArgs.supports(IN_ARG_x_dot_sg)) {
368  x_dot_sg = inArgs.get_x_dot_sg();
369  if (x_dot_sg != Teuchos::null) {
370  do_quad = true;
371  }
372  }
373  for (int i=0; i<num_p; i++) {
374  p_sg[i] = inArgs.get_p_sg(i);
375  if (p_sg[i] != Teuchos::null) {
376  do_quad = true;
377  }
378  }
379  if (outArgs.supports(OUT_ARG_f_sg)) {
380  f_sg = outArgs.get_f_sg();
381  if (f_sg != Teuchos::null)
382  f_sg->init(0.0);
383  }
384  if (outArgs.supports(OUT_ARG_W_sg)) {
385  W_sg = outArgs.get_W_sg();
386  if (W_sg != Teuchos::null)
387  W_sg->init(0.0);
388  }
389  for (int i=0; i<num_p; i++) {
390  if (!outArgs.supports(OUT_ARG_DfDp_sg, i).none()) {
391  dfdp_sg[i] = outArgs.get_DfDp_sg(i);
392  if (dfdp_sg[i].getMultiVector() != Teuchos::null)
393  dfdp_sg[i].getMultiVector()->init(0.0);
394  else if (dfdp_sg[i].getLinearOp() != Teuchos::null)
395  dfdp_sg[i].getLinearOp()->init(0.0);
396  }
397  }
398 
399  for (int i=0; i<num_g; i++) {
400  g_sg[i] = outArgs.get_g_sg(i);
401  if (g_sg[i] != Teuchos::null)
402  g_sg[i]->init(0.0);
403 
404  if (!outArgs.supports(OUT_ARG_DgDx_sg, i).none()) {
405  dgdx_sg[i] = outArgs.get_DgDx_sg(i);
406  if (dgdx_sg[i].getMultiVector() != Teuchos::null)
407  dgdx_sg[i].getMultiVector()->init(0.0);
408  else if (dgdx_sg[i].getLinearOp() != Teuchos::null)
409  dgdx_sg[i].getLinearOp()->init(0.0);
410  }
411 
412  if (!outArgs.supports(OUT_ARG_DgDx_dot_sg, i).none()) {
413  dgdx_dot_sg[i] = outArgs.get_DgDx_dot_sg(i);
414  if (dgdx_dot_sg[i].getMultiVector() != Teuchos::null)
415  dgdx_dot_sg[i].getMultiVector()->init(0.0);
416  else if (dgdx_dot_sg[i].getLinearOp() != Teuchos::null)
417  dgdx_dot_sg[i].getLinearOp()->init(0.0);
418  }
419 
420  dgdp_sg[i].resize(num_p);
421  for (int j=0; j<num_p; j++) {
422  if (!outArgs.supports(OUT_ARG_DgDp_sg, i, j).none()) {
423  dgdp_sg[i][j] = outArgs.get_DgDp_sg(i,j);
424  if (dgdp_sg[i][j].getMultiVector() != Teuchos::null)
425  dgdp_sg[i][j].getMultiVector()->init(0.0);
426  else if (dgdp_sg[i][j].getLinearOp() != Teuchos::null)
427  dgdp_sg[i][j].getLinearOp()->init(0.0);
428  }
429  }
430  }
431 
432  if (do_quad) {
433  // Get quadrature data
434  const Teuchos::Array< Teuchos::Array<double> >& quad_points =
435  quad->getQuadPoints();
436  const Teuchos::Array<double>& quad_weights =
437  quad->getQuadWeights();
438  const Teuchos::Array< Teuchos::Array<double> > & quad_values =
439  quad->getBasisAtQuadPoints();
440  const Teuchos::Array<double>& basis_norms = basis->norm_squared();
441 
442  // Perform integrations
443  for (int qp=0; qp<quad_points.size(); qp++) {
444 
445  // StieltjesPCEBasis can introduce quadrature points with zero weight
446  // Don't do those evaluations, since the model might not like the
447  // quadrature points (i.e., zero)
448  if (quad_weights[qp] == 0.0)
449  continue;
450 
451  {
452 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
453  TEUCHOS_FUNC_TIME_MONITOR_DIFF("Stokhos: SGQuadModelEvaluator -- Polynomial Evaluation",
454  PolyEvaluation);
455 #endif
456 
457  // Evaluate inputs at quadrature points
458  if (x_sg != Teuchos::null) {
459 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
460  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- X Evaluation");
461 #endif
462  x_sg->evaluate(quad_values[qp], *x_qp);
463  me_inargs.set_x(x_qp);
464  }
465  if (x_dot_sg != Teuchos::null) {
466 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
467  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- X_dot Evaluation");
468 #endif
469  x_dot_sg->evaluate(quad_values[qp], *x_dot_qp);
470  me_inargs.set_x_dot(x_qp);
471  }
472  for (int i=0; i<num_p; i++) {
473  if (p_sg[i] != Teuchos::null) {
474 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
475  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- P Evaluation");
476 #endif
477  p_sg[i]->evaluate(quad_values[qp], *(p_qp[i]));
478  me_inargs.set_p(i, p_qp[i]);
479  }
480  }
481  if (f_sg != Teuchos::null)
482  me_outargs.set_f(f_qp);
483  if (W_sg != Teuchos::null)
484  me_outargs.set_W(W_qp);
485  for (int i=0; i<num_p; i++) {
486  if (!dfdp_sg[i].isEmpty())
487  me_outargs.set_DfDp(i, dfdp_qp[i]);
488  }
489  for (int i=0; i<num_g; i++) {
490  if (g_sg[i] != Teuchos::null)
491  me_outargs.set_g(i, g_qp[i]);
492  if (!dgdx_dot_sg[i].isEmpty())
493  me_outargs.set_DgDx_dot(i, dgdx_dot_qp[i]);
494  if (!dgdx_sg[i].isEmpty())
495  me_outargs.set_DgDx(i, dgdx_qp[i]);
496  for (int j=0; j<num_p; j++)
497  if (!dgdp_sg[i][j].isEmpty())
498  me_outargs.set_DgDp(i, j, dgdp_qp[i][j]);
499  }
500 
501  }
502 
503  {
504 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
505  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- Model Evaluation");
506 #endif
507 
508  // Evaluate model at quadrature points
509  me->evalModel(me_inargs, me_outargs);
510 
511  }
512 
513  {
514 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
516  "Stokhos: SGQuadModelEvaluator -- Polynomial Integration", Integration);
517 #endif
518 
519  // Sum in results
520  if (f_sg != Teuchos::null) {
521 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
522  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- F Integration");
523 #endif
524  f_sg->sumIntoAllTerms(quad_weights[qp], quad_values[qp], basis_norms,
525  *f_qp);
526  }
527  if (W_sg != Teuchos::null) {
528 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
529  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- W Integration");
530 #endif
531  W_sg->sumIntoAllTerms(quad_weights[qp], quad_values[qp], basis_norms,
532  *W_qp);
533  }
534  for (int j=0; j<num_p; j++) {
535  if (!dfdp_sg[j].isEmpty()) {
536 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
538  "Stokhos: SGQuadModelEvaluator -- df/dp Integration");
539 #endif
540  if (dfdp_sg[j].getMultiVector() != Teuchos::null) {
541  dfdp_sg[j].getMultiVector()->sumIntoAllTerms(
542  quad_weights[qp], quad_values[qp], basis_norms,
543  *(dfdp_qp[j].getMultiVector()));
544  }
545  else if (dfdp_sg[j].getLinearOp() != Teuchos::null) {
546  dfdp_sg[j].getLinearOp()->sumIntoAllTerms(
547  quad_weights[qp], quad_values[qp], basis_norms,
548  *(dfdp_qp[j].getLinearOp()));
549  }
550  }
551  }
552  for (int i=0; i<num_g; i++) {
553  if (g_sg[i] != Teuchos::null) {
554 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
555  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- G Integration");
556 #endif
557  g_sg[i]->sumIntoAllTerms(quad_weights[qp], quad_values[qp],
558  basis_norms, *g_qp[i]);
559  }
560  if (!dgdx_dot_sg[i].isEmpty()) {
561 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
563  "Stokhos: SGQuadModelEvaluator -- dg/dx_dot Integration");
564 #endif
565  if (dgdx_dot_sg[i].getMultiVector() != Teuchos::null) {
566  dgdx_dot_sg[i].getMultiVector()->sumIntoAllTerms(
567  quad_weights[qp], quad_values[qp], basis_norms,
568  *(dgdx_dot_qp[i].getMultiVector()));
569  }
570  else if (dgdx_dot_sg[i].getLinearOp() != Teuchos::null) {
571  dgdx_dot_sg[i].getLinearOp()->sumIntoAllTerms(
572  quad_weights[qp], quad_values[qp], basis_norms,
573  *(dgdx_dot_qp[i].getLinearOp()));
574  }
575  }
576  if (!dgdx_sg[i].isEmpty()) {
577 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
579  "Stokhos: SGQuadModelEvaluator -- dg/dx Integration");
580 #endif
581  if (dgdx_sg[i].getMultiVector() != Teuchos::null) {
582  dgdx_sg[i].getMultiVector()->sumIntoAllTerms(
583  quad_weights[qp], quad_values[qp], basis_norms,
584  *(dgdx_qp[i].getMultiVector()));
585  }
586  else if (dgdx_sg[i].getLinearOp() != Teuchos::null) {
587  dgdx_sg[i].getLinearOp()->sumIntoAllTerms(
588  quad_weights[qp], quad_values[qp], basis_norms,
589  *(dgdx_qp[i].getLinearOp()));
590  }
591  }
592  for (int j=0; j<num_p; j++) {
593  if (!dgdp_sg[i][j].isEmpty()) {
594 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
596  "Stokhos: SGQuadModelEvaluator -- dg/dp Integration");
597 #endif
598  if (dgdp_sg[i][j].getMultiVector() != Teuchos::null) {
599  dgdp_sg[i][j].getMultiVector()->sumIntoAllTerms(
600  quad_weights[qp], quad_values[qp], basis_norms,
601  *(dgdp_qp[i][j].getMultiVector()));
602  }
603  else if (dgdp_sg[i][j].getLinearOp() != Teuchos::null) {
604  dgdp_sg[i][j].getLinearOp()->sumIntoAllTerms(
605  quad_weights[qp], quad_values[qp], basis_norms,
606  *(dgdp_qp[i][j].getLinearOp()));
607  }
608  }
609  }
610  }
611 
612  }
613  }
614  }
615  else {
616  // Compute the non-SG functions
617  me->evalModel(me_inargs, me_outargs);
618  }
619 }
int num_p
Number of parameter vectors.
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
int num_g
Number of response vectors.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > dgdx_dot_qp
Response derivative.
Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > dgdx_qp
Response derivative.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
InArgs createInArgs() const
Create InArgs.
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > g_qp
Response vectors.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > p_qp
Parameter vectors.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > dfdp_qp
Residual derivatives.
Teuchos::RCP< Epetra_Vector > x_dot_qp
Time derivative vector.
void resize(size_type new_size, const value_type &x=value_type())
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< Epetra_Vector > f_qp
Residual vector.
Teuchos::RCP< Epetra_Vector > x_qp
Solution vector.
Teuchos::RCP< Epetra_Operator > W_qp
W operator.
size_type size() const
Teuchos::Array< Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > > dgdp_qp
Response sensitivities.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
OutArgs createOutArgs() const
Create OutArgs.
SGQuadModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me)
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
#define TEUCHOS_FUNC_TIME_MONITOR_DIFF(FUNCNAME, DIFF)
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return observation vector map.