FEApp_SGDakotaResidualGlobalFill.cpp

Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 #include "FEApp_SGDakotaResidualGlobalFill.hpp"
00033 
00034 #if SG_ACTIVE
00035 
00036 #ifdef HAVE_DAKOTA
00037 
00038 #include "DakotaInterface.H"
00039 #include "DakotaIterator.H"
00040 #include "NonDExpansion.H"
00041 #include "DataFitSurrModel.H"
00042 
00043 FEApp::SGDakotaResidualGlobalFill::
00044 SGDakotaResidualGlobalFill(
00045       const Teuchos::RCP<const FEApp::Mesh>& elementMesh,
00046       const Teuchos::RCP<const FEApp::AbstractQuadrature>& quadRule,
00047       const Teuchos::RCP< FEApp::AbstractPDE<FEApp::SGResidualType> >& pdeEquations,
00048       const std::vector< Teuchos::RCP<FEApp::NodeBC> >& nodeBCs,
00049       bool is_transient,
00050       const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sgBasis,
00051       const Teuchos::RCP< FEApp::AbstractPDE<FEApp::ResidualType> >& resPDEEquations,
00052       const ParamVec* pvec):
00053   GlobalFill<SGResidualType>(elementMesh, quadRule, pdeEquations, nodeBCs,
00054                              is_transient),
00055   sg_basis(sgBasis),
00056   residPDE(resPDEEquations),
00057   p(pvec),
00058   sg_size(sg_basis->size()),
00059   parallel_lib(),
00060   problem_db(parallel_lib)
00061 {
00062   const char* dakota_in = "dakota_res.in";
00063   const char* dakota_out = "dakota_res.out";
00064   const char* dakota_err = "dakota_res.err";
00065   const char* dakota_restart_out = "dakota_res.restart";
00066   parallel_lib.specify_outputs_restart(dakota_out, dakota_err, NULL,
00067                                        dakota_restart_out, 0);
00068   problem_db.manage_inputs(dakota_in);
00069   problem_db.check_input();
00070   selected_strategy = Dakota::Strategy(problem_db);
00071   Dakota::Model& first_model = *(problem_db.model_list().begin());
00072   Dakota::Interface& interface = first_model.interface();
00073   elemInterface = 
00074     new FEApp::DakotaElementResidualInterface(problem_db, residPDE, *p, 
00075                                               quad, ndof);
00076   interface.assign_rep(elemInterface, false);
00077 }
00078 
00079 FEApp::SGDakotaResidualGlobalFill::
00080 ~SGDakotaResidualGlobalFill()
00081 {
00082 }
00083 
00084 void
00085 FEApp::SGDakotaResidualGlobalFill::
00086 computeGlobalFill(FEApp::AbstractInitPostOp<FEApp::SGResidualType>& initPostOp)
00087 {
00088   // Loop over elements
00089   Teuchos::RCP<const FEApp::AbstractElement> e;
00090   for (FEApp::Mesh::const_iterator eit=mesh->begin(); eit!=mesh->end(); ++eit){
00091     e = *eit;
00092 
00093     // Initialize element solution
00094     initPostOp.elementInit(*e, neqn, elem_xdot, elem_x);
00095 
00096     // Reset Dakota element interface
00097     elemInterface->reset(e, Teuchos::rcp(elem_xdot,false), 
00098                          Teuchos::rcp(&elem_x,false));
00099 
00100     // Run Dakota on this element
00101     selected_strategy.run_strategy();
00102 
00103     // Extract SG components from Dakota
00104     // We are assuming Dakota and Sacado order them the same way
00105     Dakota::Iterator& it = *(problem_db.iterator_list().begin());
00106     Dakota::NonDExpansion& nondexp = 
00107       dynamic_cast<Dakota::NonDExpansion&>(*(it.iterator_rep()));
00108     Dakota::Model& sg_model = nondexp.get_uSpaceModel();
00109     const Dakota::RealVectorArray& f_coeffs = 
00110       sg_model.approximation_coefficients();
00111     for (std::size_t i=0; i<ndof; i++) {
00112       elem_f[i].copyForWrite();
00113       elem_f[i].resize(sg_size);
00114       for (std::size_t j=0; j<sg_size; j++)
00115         elem_f[i].fastAccessCoeff(j) = f_coeffs[i][j];
00116     }
00117 
00118     // Post-process element residual
00119     initPostOp.elementPost(*e, neqn, elem_f);
00120 
00121   }
00122 
00123   // Loop over boundary conditions
00124   for (std::size_t i=0; i<bc.size(); i++) {
00125 
00126     if (bc[i]->isOwned() || bc[i]->isShared()) {
00127 
00128       // Zero out node residual
00129       for (unsigned int j=0; j<neqn; j++)
00130         node_f[j] = 0.0;
00131 
00132       // Initialize node solution
00133       initPostOp.nodeInit(*bc[i], neqn, node_xdot, node_x);
00134 
00135       // Compute node residual
00136       bc[i]->getStrategy<SGResidualType>()->evaluateResidual(node_xdot, 
00137                                                              node_x, 
00138                                                              node_f);
00139 
00140       // Post-process node residual
00141       initPostOp.nodePost(*bc[i], neqn, node_f);
00142 
00143     }
00144     
00145   }
00146 
00147   // Finalize fill
00148   initPostOp.finalizeFill();
00149 
00150 }
00151 
00152 #endif
00153 
00154 #endif

Generated on Wed May 12 21:39:33 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7