55 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_) :
71 InArgs me_inargs =
me->createInArgs();
72 num_p = me_inargs.Np();
73 if (me_inargs.supports(IN_ARG_x_dot))
75 if (me_inargs.supports(IN_ARG_x))
78 for (
int i=0; i<
num_p; i++)
82 OutArgs me_outargs =
me->createOutArgs();
83 num_g = me_outargs.Ng();
86 if (me_outargs.supports(OUT_ARG_f))
90 if (me_outargs.supports(OUT_ARG_W))
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(
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(
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));
114 for (
int i=0; i<
num_g; i++) {
121 if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_TRANS_MV_BY_ROW))
122 dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
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(
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));
136 if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_TRANS_MV_BY_ROW))
137 dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
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(
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));
153 if (me_outargs.supports(OUT_ARG_DgDp, i,
j).supports(DERIV_TRANS_MV_BY_ROW))
154 dgdp_qp[i][
j] = EpetraExt::ModelEvaluator::Derivative(
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(
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));
256 OutArgsSetup outArgs;
257 OutArgs me_outargs = me->createOutArgs();
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));
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);
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));
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));
301evalModel(
const InArgs& inArgs,
const OutArgs& outArgs)
const
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));
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));
338 bool do_quad =
false;
339 InArgs::sg_const_vector_t x_sg;
340 InArgs::sg_const_vector_t x_dot_sg;
341 Teuchos::Array<InArgs::sg_const_vector_t> p_sg(num_p);
342 OutArgs::sg_vector_t f_sg;
343 OutArgs::sg_operator_t W_sg;
344 Teuchos::Array<SGDerivative> dfdp_sg(num_p);
345 Teuchos::Array<OutArgs::sg_vector_t> g_sg(num_g);
346 Teuchos::Array<SGDerivative> dgdx_sg(num_g);
347 Teuchos::Array<SGDerivative> dgdx_dot_sg(num_g);
348 Teuchos::Array< Teuchos::Array<SGDerivative> > dgdp_sg(num_g);
349 TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_basis() == Teuchos::null,
351 "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
352 "SG basis inArg cannot be null!");
353 TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_quadrature() == Teuchos::null,
355 "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
356 "SG quadrature inArg cannot be null!");
357 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
358 inArgs.get_sg_basis();
359 Teuchos::RCP< const Stokhos::Quadrature<int,double> > quad =
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) {
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) {
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) {
379 if (outArgs.supports(OUT_ARG_f_sg)) {
380 f_sg = outArgs.get_f_sg();
381 if (f_sg != Teuchos::null)
384 if (outArgs.supports(OUT_ARG_W_sg)) {
385 W_sg = outArgs.get_W_sg();
386 if (W_sg != Teuchos::null)
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);
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)
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);
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);
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);
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();
443 for (
int qp=0; qp<quad_points.size(); qp++) {
448 if (quad_weights[qp] == 0.0)
452#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
453 TEUCHOS_FUNC_TIME_MONITOR_DIFF(
"Stokhos: SGQuadModelEvaluator -- Polynomial Evaluation",
458 if (x_sg != Teuchos::null) {
459#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
460 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: SGQuadModelEvaluator -- X Evaluation");
462 x_sg->evaluate(quad_values[qp], *x_qp);
463 me_inargs.set_x(x_qp);
465 if (x_dot_sg != Teuchos::null) {
466#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
467 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: SGQuadModelEvaluator -- X_dot Evaluation");
469 x_dot_sg->evaluate(quad_values[qp], *x_dot_qp);
470 me_inargs.set_x_dot(x_qp);
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");
477 p_sg[i]->evaluate(quad_values[qp], *(p_qp[i]));
478 me_inargs.set_p(i, p_qp[i]);
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]);
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]);
504#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
505 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: SGQuadModelEvaluator -- Model Evaluation");
509 me->evalModel(me_inargs, me_outargs);
514#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
515 TEUCHOS_FUNC_TIME_MONITOR_DIFF(
516 "Stokhos: SGQuadModelEvaluator -- Polynomial Integration", Integration);
520 if (f_sg != Teuchos::null) {
521#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
522 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: SGQuadModelEvaluator -- F Integration");
524 f_sg->sumIntoAllTerms(quad_weights[qp], quad_values[qp], basis_norms,
527 if (W_sg != Teuchos::null) {
528#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
529 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: SGQuadModelEvaluator -- W Integration");
531 W_sg->sumIntoAllTerms(quad_weights[qp], quad_values[qp], basis_norms,
534 for (
int j=0;
j<num_p;
j++) {
535 if (!dfdp_sg[
j].isEmpty()) {
536#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
537 TEUCHOS_FUNC_TIME_MONITOR(
538 "Stokhos: SGQuadModelEvaluator -- df/dp Integration");
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()));
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()));
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");
557 g_sg[i]->sumIntoAllTerms(quad_weights[qp], quad_values[qp],
558 basis_norms, *g_qp[i]);
560 if (!dgdx_dot_sg[i].isEmpty()) {
561#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
562 TEUCHOS_FUNC_TIME_MONITOR(
563 "Stokhos: SGQuadModelEvaluator -- dg/dx_dot Integration");
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()));
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()));
576 if (!dgdx_sg[i].isEmpty()) {
577#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
578 TEUCHOS_FUNC_TIME_MONITOR(
579 "Stokhos: SGQuadModelEvaluator -- dg/dx Integration");
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()));
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()));
592 for (
int j=0;
j<num_p;
j++) {
593 if (!dgdp_sg[i][
j].isEmpty()) {
594#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
595 TEUCHOS_FUNC_TIME_MONITOR(
596 "Stokhos: SGQuadModelEvaluator -- dg/dp Integration");
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()));
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()));
617 me->evalModel(me_inargs, me_outargs);