Albany: a Trilinos-based PDE code: Albany_Application.cpp Source File

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 #include "Albany_Application.hpp"
00007 #include "Albany_Utils.hpp"
00008 #include "AAdapt_AdaptationFactory.hpp"
00009 #include "Albany_ProblemFactory.hpp"
00010 #include "Albany_DiscretizationFactory.hpp"
00011 #include "Albany_ResponseFactory.hpp"
00012 #include "Epetra_LocalMap.h"
00013 #include "Stokhos_OrthogPolyBasis.hpp"
00014 #include "Teuchos_TimeMonitor.hpp"
00015 
00016 #include<string>
00017 #include "PHAL_Workset.hpp"
00018 #include "Albany_DataTypes.hpp"
00019 
00020 #include "Albany_DummyParameterAccessor.hpp"
00021 #ifdef ALBANY_CUTR
00022   #include "CUTR_CubitMeshMover.hpp"
00023   #include "STKMeshData.hpp"
00024 #endif
00025 
00026 #include "Teko_InverseFactoryOperator.hpp"
00027 #include "Teko_StridedEpetraOperator.hpp"
00028 
00029 #include "Albany_ScalarResponseFunction.hpp"
00030 
00031 #include "EpetraExt_RowMatrixOut.h"
00032 #include "EpetraExt_MultiVectorOut.h"
00033 
00034 using Teuchos::ArrayRCP;
00035 using Teuchos::RCP;
00036 using Teuchos::rcp;
00037 using Teuchos::rcp_dynamic_cast;
00038 using Teuchos::TimeMonitor;
00039 
00040 int countJac; //counter which counts instances of Jacobian (for debug output)
00041 int countRes; //counter which counts instances of residual (for debug output)
00042 
00043 
00044 Albany::Application::
00045 Application(const RCP<const Epetra_Comm>& comm_,
00046       const RCP<Teuchos::ParameterList>& params,
00047       const RCP<const Epetra_Vector>& initial_guess) :
00048   comm(comm_),
00049   out(Teuchos::VerboseObjectBase::getDefaultOStream()),
00050   physicsBasedPreconditioner(false),
00051   shapeParamsHaveBeenReset(false),
00052   morphFromInit(true), perturbBetaForDirichlets(0.0),
00053   phxGraphVisDetail(0),
00054   stateGraphVisDetail(0)
00055 {
00056   // Create parameter library
00057   paramLib = rcp(new ParamLib);
00058 
00059 #ifdef ALBANY_DEBUG
00060   int break_set = (getenv("ALBANY_BREAK") == NULL)?0:1;
00061   int env_status = 0;
00062   int length = 1;
00063   comm->SumAll(&break_set, &env_status, length);
00064   if(env_status != 0){
00065     *out << "Host and Process Ids for tasks" << std::endl;
00066     comm->Barrier();
00067     int nproc = comm->NumProc();
00068     for(int i = 0; i < nproc; i++) {
00069       if(i == comm->MyPID()) {
00070         char buf[80];
00071         char hostname[80]; gethostname(hostname, sizeof(hostname));
00072         sprintf(buf, "Host: %s   PID: %d", hostname, getpid());
00073         *out << buf << std::endl;
00074         std::cout.flush();
00075         sleep(1);
00076       }
00077       comm->Barrier();
00078     }
00079     if(comm->MyPID() == 0) {
00080       char go = ' ';
00081       std::cout << "\n";
00082       std::cout << "** Client has paused because the environment variable ALBANY_BREAK has been set.\n";
00083       std::cout << "** You may attach a debugger to processes now.\n";
00084       std::cout << "**\n";
00085       std::cout << "** Enter a character (not whitespace), then <Return> to continue. > "; std::cout.flush();
00086       std::cin >> go;
00087       std::cout << "\n** Now pausing for 3 seconds.\n"; std::cout.flush();
00088     }
00089     sleep(3);
00090   }
00091   comm->Barrier();
00092 #endif
00093 
00094   // Create problem object
00095   RCP<Teuchos::ParameterList> problemParams =
00096     Teuchos::sublist(params, "Problem", true);
00097   Albany::ProblemFactory problemFactory(problemParams, paramLib, comm);
00098   problem = problemFactory.create();
00099 
00100   // Validate Problem parameters against list for this specific problem
00101   problemParams->validateParameters(*(problem->getValidProblemParameters()),0);
00102 
00103   // Save the solution method to be used
00104   std::string solutionMethod = problemParams->get("Solution Method", "Steady");
00105   if(solutionMethod == "Steady")
00106     solMethod = Steady;
00107   else if(solutionMethod == "Continuation")
00108     solMethod = Continuation;
00109   else if(solutionMethod == "Transient")
00110     solMethod = Transient;
00111   else if(solutionMethod == "Eigensolve")
00112     solMethod = Eigensolve;
00113   else
00114     TEUCHOS_TEST_FOR_EXCEPTION(true,
00115             std::logic_error, "Solution Method must be Steady, Transient, "
00116             << "Continuation, or Eigensolve not : " << solutionMethod);
00117 
00118   // Register shape parameters for manipulation by continuation/optimization
00119   if (problemParams->get("Enable Cubit Shape Parameters",false)) {
00120 #ifdef ALBANY_CUTR
00121     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
00122     meshMover = rcp(new CUTR::CubitMeshMover
00123           (problemParams->get<std::string>("Cubit Base Filename")));
00124 
00125     meshMover->getShapeParams(shapeParamNames, shapeParams);
00126     *out << "SSS : Registering " << shapeParams.size() << " Shape Parameters" << std::endl;
00127 
00128     registerShapeParameters();
00129 
00130 #else
00131   TEUCHOS_TEST_FOR_EXCEPTION(problemParams->get("Enable Cubit Shape Parameters",false), std::logic_error,
00132            "Cubit requested but not Compiled in!");
00133 #endif
00134   }
00135 
00136   determinePiroSolver(params);
00137 
00138   physicsBasedPreconditioner = problemParams->get("Use Physics-Based Preconditioner",false);
00139   if (physicsBasedPreconditioner)
00140     tekoParams = Teuchos::sublist(problemParams, "Teko", true);
00141 
00142   // Create debug output object
00143   RCP<Teuchos::ParameterList> debugParams =
00144     Teuchos::sublist(params, "Debug Output", true);
00145   writeToMatrixMarketJac = debugParams->get("Write Jacobian to MatrixMarket", 0);
00146   writeToMatrixMarketRes = debugParams->get("Write Residual to MatrixMarket", 0);
00147   writeToCoutJac = debugParams->get("Write Jacobian to Standard Output", 0);
00148   writeToCoutRes = debugParams->get("Write Residual to Standard Output", 0);
00149   //the above 4 parameters cannot have values < -1
00150   if (writeToMatrixMarketJac < -1)  {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00151           std::endl << "Error in Albany::Application constructor:  " <<
00152           "Invalid Parameter Write Jacobian to MatrixMarket.  Acceptable values are -1, 0, 1, 2, ... " << std::endl);}
00153   if (writeToMatrixMarketRes < -1)  {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00154           std::endl << "Error in Albany::Application constructor:  " <<
00155           "Invalid Parameter Write Residual to MatrixMarket.  Acceptable values are -1, 0, 1, 2, ... " << std::endl);}
00156   if (writeToCoutJac < -1)  {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00157           std::endl << "Error in Albany::Application constructor:  " <<
00158           "Invalid Parameter Write Jacobian to Standard Output.  Acceptable values are -1, 0, 1, 2, ... " << std::endl);}
00159   if (writeToCoutRes < -1)  {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00160           std::endl << "Error in Albany::Application constructor:  " <<
00161           "Invalid Parameter Write Residual to Standard Output.  Acceptable values are -1, 0, 1, 2, ... " << std::endl);}
00162   if (writeToMatrixMarketJac != 0 || writeToCoutJac != 0 )
00163      countJac = 0; //initiate counter that counts instances of Jacobian matrix to 0
00164   if (writeToMatrixMarketRes != 0 || writeToCoutRes != 0)
00165      countRes = 0; //initiate counter that counts instances of Jacobian matrix to 0
00166 
00167   // Create discretization object
00168 
00169   Albany::DiscretizationFactory discFactory(params, comm);
00170 
00171 #ifdef ALBANY_CUTR
00172   discFactory.setMeshMover(meshMover);
00173 #endif
00174 
00175   // Get mesh specification object: worksetSize, cell topology, etc
00176   ArrayRCP<RCP<Albany::MeshSpecsStruct> > meshSpecs =
00177     discFactory.createMeshSpecs();
00178 
00179   problem->buildProblem(meshSpecs, stateMgr);
00180 
00181   // Construct responses
00182   // This really needs to happen after the discretization is created for
00183   // distributed responses, but currently it can't be moved because there
00184   // are responses that setup states, which has to happen before the
00185   // discreatization is created.  We will delay setup of the distributed
00186   // responses to deal with this temporarily.
00187   Teuchos::ParameterList& responseList =
00188     problemParams->sublist("Response Functions");
00189   ResponseFactory responseFactory(Teuchos::rcp(this,false), problem, meshSpecs,
00190           Teuchos::rcp(&stateMgr,false));
00191   responses = responseFactory.createResponseFunctions(responseList);
00192 
00193   // Build state field manager
00194   sfm.resize(meshSpecs.size());
00195   Teuchos::RCP<PHX::DataLayout> dummy =
00196     Teuchos::rcp(new PHX::MDALayout<Dummy>(0));
00197   for (int ps=0; ps<meshSpecs.size(); ps++) {
00198     std::string elementBlockName = meshSpecs[ps]->ebName;
00199     std::vector<std::string>responseIDs_to_require =
00200       stateMgr.getResidResponseIDsToRequire(elementBlockName);
00201     sfm[ps] = Teuchos::rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
00202     Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> > tags =
00203       problem->buildEvaluators(*sfm[ps], *meshSpecs[ps], stateMgr,
00204              BUILD_STATE_FM, Teuchos::null);
00205     std::vector<std::string>::const_iterator it;
00206     for (it = responseIDs_to_require.begin();
00207    it != responseIDs_to_require.end();
00208    it++) {
00209       const std::string& responseID = *it;
00210       PHX::Tag<PHAL::AlbanyTraits::Residual::ScalarT> res_response_tag(
00211   responseID, dummy);
00212       sfm[ps]->requireField<PHAL::AlbanyTraits::Residual>(res_response_tag);
00213     }
00214     sfm[ps]->postRegistrationSetup("");
00215   }
00216 
00217   // Create the full mesh
00218   neq = problem->numEquations();
00219   disc = discFactory.createDiscretization(neq, stateMgr.getStateInfoStruct(), 
00220                                           problem->getFieldRequirements(),
00221                                           problem->getNullSpace());
00222 
00223   solMgr = rcp(new AAdapt::AdaptiveSolutionManager(params, disc, initial_guess));
00224 
00225   // Now that space is allocated in STK for state fields, initialize states
00226   stateMgr.setStateArrays(disc);
00227 
00228   // Now setup response functions (see note above)
00229   for (int i=0; i<responses.size(); i++)
00230     responses[i]->setup();
00231 
00232   // Set up memory for workset
00233   fm = problem->getFieldManager();
00234   TEUCHOS_TEST_FOR_EXCEPTION(fm==Teuchos::null, std::logic_error,
00235            "getFieldManager not implemented!!!");
00236   dfm = problem->getDirichletFieldManager();
00237   nfm = problem->getNeumannFieldManager();
00238 
00239   if (comm->MyPID()==0) {
00240     phxGraphVisDetail= problemParams->get("Phalanx Graph Visualization Detail", 0);
00241     stateGraphVisDetail= phxGraphVisDetail;
00242   }
00243 
00244   *out << "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n"
00245        << " Sacado ParameterLibrary has been initialized:\n "
00246        << *paramLib
00247        << "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n"
00248        << std::endl;
00249 
00250   ignore_residual_in_jacobian =
00251     problemParams->get("Ignore Residual In Jacobian", false);
00252 
00253   perturbBetaForDirichlets = problemParams->get("Perturb Dirichlet",0.0);
00254 
00255   is_adjoint =
00256     problemParams->get("Solve Adjoint", false);
00257 
00258   // For backward compatibility, use any value at the old location of the "Compute Sensitivity" flag
00259   // as a default value for the new flag location when the latter has been left undefined
00260   const std::string sensitivityToken = "Compute Sensitivities";
00261   const Teuchos::Ptr<const bool> oldSensitivityFlag(problemParams->getPtr<bool>(sensitivityToken));
00262   if (Teuchos::nonnull(oldSensitivityFlag)) {
00263     Teuchos::ParameterList &solveParams = params->sublist("Piro").sublist("Analysis").sublist("Solve");
00264     solveParams.get(sensitivityToken, *oldSensitivityFlag);
00265   }
00266 
00267 #ifdef ALBANY_MOR
00268   if(disc->supportsMOR())
00269     morFacade = createMORFacade(disc, problemParams);
00270 #endif
00271 
00272 /*
00273  * Initialize mesh adaptation features
00274  */
00275 
00276   if(solMgr->hasAdaptation()){
00277 
00278     solMgr->buildAdaptiveProblem(paramLib, stateMgr, comm);
00279 
00280   }
00281 
00282 }
00283 
00284 Albany::Application::
00285 ~Application()
00286 {
00287 #ifdef ALBANY_DEBUG
00288   *out << "Calling destructor for Albany_Application" << std::endl;
00289 #endif
00290 }
00291 
00292 RCP<Albany::AbstractDiscretization>
00293 Albany::Application::
00294 getDiscretization() const
00295 {
00296   return disc;
00297 }
00298 
00299 RCP<Albany::AbstractProblem>
00300 Albany::Application::
00301 getProblem() const
00302 {
00303   return problem;
00304 }
00305 
00306 RCP<const Epetra_Comm>
00307 Albany::Application::
00308 getComm() const
00309 {
00310   return comm;
00311 }
00312 
00313 RCP<const Epetra_Map>
00314 Albany::Application::
00315 getMap() const
00316 {
00317   return disc->getMap();
00318 }
00319 
00320 RCP<const Epetra_CrsGraph>
00321 Albany::Application::
00322 getJacobianGraph() const
00323 {
00324   return disc->getJacobianGraph();
00325 }
00326 
00327 
00328 RCP<Epetra_Operator>
00329 Albany::Application::
00330 getPreconditioner()
00331 {
00332    //inverseLib = Teko::InverseLibrary::buildFromStratimikos();
00333    inverseLib = Teko::InverseLibrary::buildFromParameterList(tekoParams->sublist("Inverse Factory Library"));
00334    inverseLib->PrintAvailableInverses(*out);
00335 
00336    inverseFac = inverseLib->getInverseFactory(tekoParams->get("Preconditioner Name","Amesos"));
00337 
00338    // get desired blocking of unknowns
00339    std::stringstream ss;
00340    ss << tekoParams->get<std::string>("Unknown Blocking");
00341 
00342    // figure out the decomposition requested by the string
00343    unsigned int num=0,sum=0;
00344    while(not ss.eof()) {
00345       ss >> num;
00346       TEUCHOS_ASSERT(num>0);
00347       sum += num;
00348       blockDecomp.push_back(num);
00349    }
00350    TEUCHOS_ASSERT(neq==sum);
00351 
00352    return rcp(new Teko::Epetra::InverseFactoryOperator(inverseFac));
00353 }
00354 
00355 RCP<ParamLib>
00356 Albany::Application::
00357 getParamLib()
00358 {
00359   return paramLib;
00360 }
00361 
00362 int
00363 Albany::Application::
00364 getNumResponses() const {
00365   return responses.size();
00366 }
00367 
00368 Teuchos::RCP<Albany::AbstractResponseFunction>
00369 Albany::Application::
00370 getResponse(int i) const
00371 {
00372   return responses[i];
00373 }
00374 
00375 bool
00376 Albany::Application::
00377 suppliesPreconditioner() const
00378 {
00379   return physicsBasedPreconditioner;
00380 }
00381 
00382 RCP<Stokhos::OrthogPolyExpansion<int,double> >
00383 Albany::Application::
00384 getStochasticExpansion()
00385 {
00386   return sg_expansion;
00387 }
00388 
00389 #ifdef ALBANY_SG_MP
00390 void
00391 Albany::Application::
00392 init_sg(const RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis,
00393   const RCP<const Stokhos::Quadrature<int,double> >& quad,
00394   const RCP<Stokhos::OrthogPolyExpansion<int,double> >& expansion,
00395   const RCP<const EpetraExt::MultiComm>& multiComm)
00396 {
00397 
00398   // Setup stohastic Galerkin
00399   sg_basis = basis;
00400   sg_quad = quad;
00401   sg_expansion = expansion;
00402   product_comm = multiComm;
00403 
00404   if (sg_overlapped_x == Teuchos::null) {
00405     sg_overlap_map =
00406       rcp(new Epetra_LocalMap(sg_basis->size(), 0,
00407             product_comm->TimeDomainComm()));
00408     sg_overlapped_x =
00409       rcp(new Stokhos::EpetraVectorOrthogPoly(
00410       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
00411     sg_overlapped_xdot =
00412   rcp(new Stokhos::EpetraVectorOrthogPoly(
00413         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
00414     sg_overlapped_xdotdot =
00415   rcp(new Stokhos::EpetraVectorOrthogPoly(
00416         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
00417     sg_overlapped_f =
00418       rcp(new Stokhos::EpetraVectorOrthogPoly(
00419       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
00420     // Delay creation of sg_overlapped_jac until needed
00421   }
00422 
00423   // Initialize responses
00424   for (int i=0; i<responses.size(); i++)
00425     responses[i]->init_sg(basis, quad, expansion, multiComm);
00426 }
00427 #endif //ALBANY_SG_MP
00428 
00429 void
00430 Albany::Application::
00431 computeGlobalResidual(const double current_time,
00432           const Epetra_Vector* xdot,
00433           const Epetra_Vector* xdotdot,
00434           const Epetra_Vector& x,
00435           const Teuchos::Array<ParamVec>& p,
00436           Epetra_Vector& f)
00437 {
00438   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Residual");
00439 
00440   postRegSetup("Residual");
00441 
00442   // Load connectivity map and coordinates
00443   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
00444         wsElNodeEqID = disc->getWsElNodeEqID();
00445   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
00446         coords = disc->getCoords();
00447   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
00448   //      sHeight = disc->getSurfaceHeight();
00449   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
00450   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
00451 
00452   int numWorksets = wsElNodeEqID.size();
00453 
00454   Teuchos::RCP<Epetra_Vector>& overlapped_f = solMgr->get_overlapped_f();
00455   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
00456 
00457   // Scatter x and xdot to the overlapped distrbution
00458   solMgr->scatterX(x, xdot, xdotdot);
00459 
00460   // Set parameters
00461   for (int i=0; i<p.size(); i++)
00462     for (unsigned int j=0; j<p[i].size(); j++)
00463       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
00464 
00465   // Mesh motion needs to occur here on the global mesh befor
00466   // it is potentially carved into worksets.
00467 #ifdef ALBANY_CUTR
00468   static int first=true;
00469   if (shapeParamsHaveBeenReset) {
00470     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
00471 
00472 *out << " Calling moveMesh with params: " << std::setprecision(8);
00473  for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
00474 *out << std::endl;
00475     meshMover->moveMesh(shapeParams, morphFromInit);
00476     coords = disc->getCoords();
00477     shapeParamsHaveBeenReset = false;
00478   }
00479 #endif
00480 
00481 
00482   // Zero out overlapped residual
00483   overlapped_f->PutScalar(0.0);
00484   f.PutScalar(0.0);
00485 
00486   // Set data in Workset struct, and perform fill via field manager
00487   {
00488     PHAL::Workset workset;
00489 
00490     if (!paramLib->isParameter("Time"))
00491 //      loadBasicWorksetInfo( workset, overlapped_x, overlapped_xdot, current_time );
00492       loadBasicWorksetInfo( workset, current_time );
00493     else
00494 //      loadBasicWorksetInfo( workset, overlapped_x, overlapped_xdot,
00495       loadBasicWorksetInfo( workset,
00496           paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time") );
00497 
00498     workset.f        = overlapped_f;
00499 
00500     for (int ws=0; ws < numWorksets; ws++) {
00501 
00502       loadWorksetBucketInfo<PHAL::AlbanyTraits::Residual>(workset, ws);
00503 
00504       // FillType template argument used to specialize Sacado
00505       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Residual>(workset);
00506       if (nfm!=Teuchos::null)
00507          nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Residual>(workset);
00508     }
00509   }
00510 
00511   f.Export(*overlapped_f, *exporter, Add);
00512 
00513   disc->setResidualField(f);
00514 
00515   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
00516   if (dfm!=Teuchos::null) {
00517     PHAL::Workset workset;
00518 
00519     workset.f = Teuchos::rcpFromRef(f);
00520     loadWorksetNodesetInfo(workset);
00521     workset.x = Teuchos::rcpFromRef(x);
00522     if ( paramLib->isParameter("Time") )
00523       workset.current_time = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
00524     else
00525       workset.current_time = current_time;
00526     if (xdot != NULL) workset.transientTerms = true;
00527     if (xdotdot != NULL) workset.accelerationTerms = true;
00528 
00529     // FillType template argument used to specialize Sacado
00530     dfm->evaluateFields<PHAL::AlbanyTraits::Residual>(workset);
00531   }
00532   //cout << "Global Resid f\n" << f << std::endl;
00533   //std::cout << "Global Soln x\n" << x << std::endl;
00534 
00535   //Debut output
00536   if (writeToMatrixMarketRes != 0) { //If requesting writing to MatrixMarket of residual...
00537     char name[100];  //create string for file name
00538     if (writeToMatrixMarketRes == -1) { //write residual to MatrixMarket every time it arises
00539        sprintf(name, "rhs%i.mm", countRes);
00540        EpetraExt::MultiVectorToMatrixMarketFile(name, f);
00541     }
00542     else {
00543       if (countRes == writeToMatrixMarketRes) { //write residual only at requested count#
00544         sprintf(name, "rhs%i.mm", countRes);
00545         EpetraExt::MultiVectorToMatrixMarketFile(name, f);
00546       }
00547     }
00548   }
00549   if (writeToCoutRes != 0) { //If requesting writing of residual to cout...
00550     if (writeToCoutRes == -1) { //cout residual time it arises
00551        std::cout << "Global Residual #" << countRes << ": " << std::endl;
00552        std::cout << f << std::endl;
00553     }
00554     else {
00555       if (countRes == writeToCoutRes) { //cout residual only at requested count#
00556         std::cout << "Global Residual #" << countRes << ": " << std::endl;
00557         std::cout << f << std::endl;
00558       }
00559     }
00560   }
00561   if (writeToMatrixMarketRes != 0 || writeToCoutRes != 0)
00562     countRes++;  //increment residual counter
00563 }
00564 
00565 void
00566 Albany::Application::
00567 computeGlobalJacobian(const double alpha,
00568           const double beta,
00569           const double omega,
00570           const double current_time,
00571           const Epetra_Vector* xdot,
00572           const Epetra_Vector* xdotdot,
00573           const Epetra_Vector& x,
00574           const Teuchos::Array<ParamVec>& p,
00575           Epetra_Vector* f,
00576           Epetra_CrsMatrix& jac)
00577 {
00578   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Jacobian");
00579 
00580   postRegSetup("Jacobian");
00581 
00582   // Load connectivity map and coordinates
00583   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
00584         wsElNodeEqID = disc->getWsElNodeEqID();
00585   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
00586         coords = disc->getCoords();
00587   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
00588   //      sHeight = disc->getSurfaceHeight();
00589   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
00590   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
00591 
00592   int numWorksets = wsElNodeEqID.size();
00593 
00594   Teuchos::RCP<Epetra_Vector>& overlapped_f = solMgr->get_overlapped_f();
00595   Teuchos::RCP<Epetra_CrsMatrix>& overlapped_jac = solMgr->get_overlapped_jac();
00596   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
00597 
00598   // Scatter x and xdot to the overlapped distribution
00599   solMgr->scatterX(x, xdot, xdotdot);
00600 
00601   // Set parameters
00602   for (int i=0; i<p.size(); i++)
00603     for (unsigned int j=0; j<p[i].size(); j++)
00604       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
00605 
00606 #ifdef ALBANY_CUTR
00607   if (shapeParamsHaveBeenReset) {
00608     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
00609 
00610 *out << " Calling moveMesh with params: " << std::setprecision(8);
00611  for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
00612 *out << std::endl;
00613     meshMover->moveMesh(shapeParams, morphFromInit);
00614     coords = disc->getCoords();
00615     shapeParamsHaveBeenReset = false;
00616   }
00617 #endif
00618 
00619   // Zero out overlapped residual
00620   if (f != NULL) {
00621     overlapped_f->PutScalar(0.0);
00622     f->PutScalar(0.0);
00623   }
00624 
00625   // Zero out Jacobian
00626   overlapped_jac->PutScalar(0.0);
00627   jac.PutScalar(0.0);
00628 
00629   // Set data in Workset struct, and perform fill via field manager
00630   {
00631     PHAL::Workset workset;
00632     if (!paramLib->isParameter("Time"))
00633 //      loadBasicWorksetInfo( workset, overlapped_x, overlapped_xdot, current_time );
00634       loadBasicWorksetInfo( workset, current_time );
00635     else
00636 //      loadBasicWorksetInfo( workset, overlapped_x, overlapped_xdot,
00637       loadBasicWorksetInfo( workset,
00638           paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time") );
00639 
00640     workset.f        = overlapped_f;
00641     workset.Jac      = overlapped_jac;
00642     loadWorksetJacobianInfo(workset, alpha, beta, omega);
00643 
00644     for (int ws=0; ws < numWorksets; ws++) {
00645       loadWorksetBucketInfo<PHAL::AlbanyTraits::Jacobian>(workset, ws);
00646 
00647       // FillType template argument used to specialize Sacado
00648       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Jacobian>(workset);
00649       if (nfm!=Teuchos::null)
00650         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Jacobian>(workset);
00651     }
00652   }
00653 
00654   { TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Jacobian Export");
00655   // Assemble global residual
00656   if (f != NULL)
00657     f->Export(*overlapped_f, *exporter, Add);
00658 
00659   // Assemble global Jacobian
00660   jac.Export(*overlapped_jac, *exporter, Add);
00661   } // End timer
00662 
00663   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
00664   if (dfm!=Teuchos::null) {
00665     PHAL::Workset workset;
00666 
00667     workset.f = rcp(f,false);
00668     workset.Jac = Teuchos::rcpFromRef(jac);
00669     workset.m_coeff = alpha;
00670     workset.n_coeff = omega;
00671     workset.j_coeff = beta;
00672 
00673     if ( paramLib->isParameter("Time") )
00674       workset.current_time = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
00675     else
00676       workset.current_time = current_time;
00677 
00678     if (beta==0.0 && perturbBetaForDirichlets>0.0) workset.j_coeff = perturbBetaForDirichlets;
00679 
00680     workset.x = Teuchos::rcpFromRef(x);
00681     if (xdot != NULL) workset.transientTerms = true;
00682     if (xdotdot != NULL) workset.accelerationTerms = true;
00683 
00684     loadWorksetNodesetInfo(workset);
00685 
00686     // FillType template argument used to specialize Sacado
00687     dfm->evaluateFields<PHAL::AlbanyTraits::Jacobian>(workset);
00688   }
00689 
00690   jac.FillComplete(true);
00691   //std::cout << "f " << *f << std::endl;;
00692   //std::cout << "J " << jac << std::endl;;
00693 
00694   //Debut output
00695   if (writeToMatrixMarketJac != 0) { //If requesting writing to MatrixMarket of Jacobian...
00696     char name[100];  //create string for file name
00697     if (writeToMatrixMarketJac == -1) { //write jacobian to MatrixMarket every time it arises
00698        sprintf(name, "jac%i.mm", countJac);
00699        EpetraExt::RowMatrixToMatrixMarketFile(name, jac);
00700     }
00701     else {
00702       if (countJac == writeToMatrixMarketJac) { //write jacobian only at requested count#
00703         sprintf(name, "jac%i.mm", countJac);
00704         EpetraExt::RowMatrixToMatrixMarketFile(name, jac);
00705       }
00706     }
00707   }
00708   if (writeToCoutJac != 0) { //If requesting writing Jacobian to standard output (cout)...
00709     if (writeToCoutJac == -1) { //cout jacobian every time it arises
00710        std::cout << "Global Jacobian #" << countJac << ": " << std::endl;
00711        std::cout << jac << std::endl;
00712     }
00713     else {
00714       if (countJac == writeToCoutJac) { //cout jacobian only at requested count#
00715        std::cout << "Global Jacobian #" << countJac << ": " << std::endl;
00716        std::cout << jac << std::endl;
00717       }
00718     }
00719   }
00720   if (writeToMatrixMarketJac != 0 || writeToCoutJac != 0)
00721     countJac++; //increment Jacobian counter
00722 
00723 }
00724 
00725 void
00726 Albany::Application::
00727 computeGlobalPreconditioner(const RCP<Epetra_CrsMatrix>& jac,
00728           const RCP<Epetra_Operator>& prec)
00729 {
00730   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Precond");
00731 
00732   *out << "Computing WPrec by Teko" << std::endl;
00733 
00734   RCP<Teko::Epetra::InverseFactoryOperator> blockPrec
00735     = rcp_dynamic_cast<Teko::Epetra::InverseFactoryOperator>(prec);
00736 
00737   blockPrec->initInverse();
00738 
00739   wrappedJac = buildWrappedOperator(jac, wrappedJac);
00740   blockPrec->rebuildInverseOperator(wrappedJac);
00741 }
00742 
00743 void
00744 Albany::Application::
00745 computeGlobalTangent(const double alpha,
00746          const double beta,
00747          const double omega,
00748          const double current_time,
00749          bool sum_derivs,
00750          const Epetra_Vector* xdot,
00751          const Epetra_Vector* xdotdot,
00752          const Epetra_Vector& x,
00753          const Teuchos::Array<ParamVec>& par,
00754          ParamVec* deriv_par,
00755          const Epetra_MultiVector* Vx,
00756          const Epetra_MultiVector* Vxdot,
00757          const Epetra_MultiVector* Vxdotdot,
00758          const Epetra_MultiVector* Vp,
00759          Epetra_Vector* f,
00760          Epetra_MultiVector* JV,
00761          Epetra_MultiVector* fp)
00762 {
00763   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Tangent");
00764 
00765   postRegSetup("Tangent");
00766 
00767   // Load connectivity map and coordinates
00768   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
00769         wsElNodeEqID = disc->getWsElNodeEqID();
00770   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
00771         coords = disc->getCoords();
00772   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
00773   //       sHeight = disc->getSurfaceHeight();
00774   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
00775   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
00776 
00777   int numWorksets = wsElNodeEqID.size();
00778 
00779   Teuchos::RCP<Epetra_Vector>& overlapped_f = solMgr->get_overlapped_f();
00780   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
00781   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
00782 
00783   // Scatter x and xdot to the overlapped distribution
00784   solMgr->scatterX(x, xdot, xdotdot);
00785 
00786   // Scatter Vx dot the overlapped distribution
00787   RCP<Epetra_MultiVector> overlapped_Vx;
00788   if (Vx != NULL) {
00789     overlapped_Vx =
00790       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
00791             Vx->NumVectors()));
00792     overlapped_Vx->Import(*Vx, *importer, Insert);
00793   }
00794 
00795   // Scatter Vx dot the overlapped distribution
00796   RCP<Epetra_MultiVector> overlapped_Vxdot;
00797   if (Vxdot != NULL) {
00798     overlapped_Vxdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdot->NumVectors()));
00799     overlapped_Vxdot->Import(*Vxdot, *importer, Insert);
00800   }
00801   RCP<Epetra_MultiVector> overlapped_Vxdotdot;
00802   if (Vxdotdot != NULL) {
00803     overlapped_Vxdotdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdotdot->NumVectors()));
00804     overlapped_Vxdotdot->Import(*Vxdotdot, *importer, Insert);
00805   }
00806 
00807   // Set parameters
00808   for (int i=0; i<par.size(); i++)
00809     for (unsigned int j=0; j<par[i].size(); j++)
00810       par[i][j].family->setRealValueForAllTypes(par[i][j].baseValue);
00811 
00812   RCP<const Epetra_MultiVector > vp = rcp(Vp, false);
00813   RCP<ParamVec> params = rcp(deriv_par, false);
00814 
00815   // Zero out overlapped residual
00816   RCP<Epetra_Vector> overlapped_ff;
00817   if (f != NULL) {
00818     overlapped_ff = overlapped_f;
00819     overlapped_ff->PutScalar(0.0);
00820     f->PutScalar(0.0);
00821   }
00822 
00823   RCP<Epetra_MultiVector> overlapped_JV;
00824   if (JV != NULL) {
00825     overlapped_JV =
00826       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
00827             JV->NumVectors()));
00828     overlapped_JV->PutScalar(0.0);
00829     JV->PutScalar(0.0);
00830   }
00831 
00832   RCP<Epetra_MultiVector> overlapped_fp;
00833   if (fp != NULL) {
00834     overlapped_fp =
00835       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
00836             fp->NumVectors()));
00837     overlapped_fp->PutScalar(0.0);
00838     fp->PutScalar(0.0);
00839   }
00840 
00841   // Number of x & xdot tangent directions
00842   int num_cols_x = 0;
00843   if (Vx != NULL)
00844     num_cols_x = Vx->NumVectors();
00845   else if (Vxdot != NULL)
00846     num_cols_x = Vxdot->NumVectors();
00847   else if (Vxdotdot != NULL)
00848     num_cols_x = Vxdotdot->NumVectors();
00849 
00850   // Number of parameter tangent directions
00851   int num_cols_p = 0;
00852   if (params != Teuchos::null) {
00853     if (Vp != NULL)
00854       num_cols_p = Vp->NumVectors();
00855     else
00856       num_cols_p = params->size();
00857   }
00858 
00859   // Whether x and param tangent components are added or separate
00860   int param_offset = 0;
00861   if (!sum_derivs)
00862     param_offset = num_cols_x;  // offset of parameter derivs in deriv array
00863 
00864   TEUCHOS_TEST_FOR_EXCEPTION(sum_derivs &&
00865            (num_cols_x != 0) &&
00866            (num_cols_p != 0) &&
00867            (num_cols_x != num_cols_p),
00868            std::logic_error,
00869            "Seed matrices Vx and Vp must have the same number " <<
00870            " of columns when sum_derivs is true and both are "
00871            << "non-null!" << std::endl);
00872 
00873   // Initialize
00874   if (params != Teuchos::null) {
00875     TanFadType p;
00876     int num_cols_tot = param_offset + num_cols_p;
00877     for (unsigned int i=0; i<params->size(); i++) {
00878       p = TanFadType(num_cols_tot, (*params)[i].baseValue);
00879       if (Vp != NULL)
00880         for (int k=0; k<num_cols_p; k++)
00881           p.fastAccessDx(param_offset+k) = (*Vp)[k][i];
00882       else
00883         p.fastAccessDx(param_offset+i) = 1.0;
00884       (*params)[i].family->setValue<PHAL::AlbanyTraits::Tangent>(p);
00885     }
00886   }
00887 
00888   // Begin shape optimization logic
00889   ArrayRCP<ArrayRCP<double> > coord_derivs;
00890   // ws, sp, cell, node, dim
00891   ArrayRCP<ArrayRCP<ArrayRCP<ArrayRCP<ArrayRCP<double> > > > > ws_coord_derivs;
00892   ws_coord_derivs.resize(coords.size());
00893   std::vector<int> coord_deriv_indices;
00894 #ifdef ALBANY_CUTR
00895   if (shapeParamsHaveBeenReset) {
00896     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
00897 
00898      int num_sp = 0;
00899      std::vector<int> shape_param_indices;
00900 
00901      // Find any shape params from param list
00902      for (unsigned int i=0; i<params->size(); i++) {
00903        for (unsigned int j=0; j<shapeParamNames.size(); j++) {
00904          if ((*params)[i].family->getName() == shapeParamNames[j]) {
00905            num_sp++;
00906            coord_deriv_indices.resize(num_sp);
00907            shape_param_indices.resize(num_sp);
00908            coord_deriv_indices[num_sp-1] = i;
00909            shape_param_indices[num_sp-1] = j;
00910          }
00911        }
00912      }
00913 
00914     TEUCHOS_TEST_FOR_EXCEPTION( Vp != NULL, std::logic_error,
00915         "Derivatives with respect to a vector of shape\n " <<
00916         "parameters has not been implemented. Need to write\n" <<
00917         "directional derivative perturbation through meshMover!" <<
00918         std::endl);
00919 
00920      // Compute FD derivs of coordinate vector w.r.t. shape params
00921      double eps = 1.0e-4;
00922      double pert;
00923      coord_derivs.resize(num_sp);
00924      for (int ws=0; ws<coords.size(); ws++)  ws_coord_derivs[ws].resize(num_sp);
00925      for (int i=0; i<num_sp; i++) {
00926 *out << "XXX perturbing parameter " << coord_deriv_indices[i]
00927      << " which is shapeParam # " << shape_param_indices[i]
00928      << " with name " <<  shapeParamNames[shape_param_indices[i]]
00929      << " which should equal " << (*params)[coord_deriv_indices[i]].family->getName() << std::endl;
00930 
00931      pert = (fabs(shapeParams[shape_param_indices[i]]) + 1.0e-2) * eps;
00932 
00933        shapeParams[shape_param_indices[i]] += pert;
00934 *out << " Calling moveMesh with params: " << std::setprecision(8);
00935 for (unsigned int ii=0; ii<shapeParams.size(); ii++) *out << shapeParams[ii] << "  ";
00936 *out << std::endl;
00937        meshMover->moveMesh(shapeParams, morphFromInit);
00938        for (int ws=0; ws<coords.size(); ws++) {  //worset
00939          ws_coord_derivs[ws][i].resize(coords[ws].size());
00940          for (int e=0; e<coords[ws].size(); e++) { //cell
00941            ws_coord_derivs[ws][i][e].resize(coords[ws][e].size());
00942            for (int j=0; j<coords[ws][e].size(); j++) { //node
00943              ws_coord_derivs[ws][i][e][j].resize(disc->getNumDim());
00944              for (int d=0; d<disc->getNumDim(); d++)  //node
00945                 ws_coord_derivs[ws][i][e][j][d] = coords[ws][e][j][d];
00946        } } } }
00947 
00948        shapeParams[shape_param_indices[i]] -= pert;
00949      }
00950 *out << " Calling moveMesh with params: " << std::setprecision(8);
00951 for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
00952 *out << std::endl;
00953      meshMover->moveMesh(shapeParams, morphFromInit);
00954      coords = disc->getCoords();
00955 
00956      for (int i=0; i<num_sp; i++) {
00957        for (int ws=0; ws<coords.size(); ws++)  //worset
00958          for (int e=0; e<coords[ws].size(); e++)  //cell
00959            for (int j=0; j<coords[ws][i].size(); j++)  //node
00960              for (int d=0; d<disc->getNumDim; d++)  //node
00961                 ws_coord_derivs[ws][i][e][j][d] = (ws_coord_derivs[ws][i][e][j][d] - coords[ws][e][j][d]) / pert;
00962        }
00963      }
00964      shapeParamsHaveBeenReset = false;
00965   }
00966   // End shape optimization logic
00967 #endif
00968 
00969   // Set data in Workset struct, and perform fill via field manager
00970   {
00971     PHAL::Workset workset;
00972     if (!paramLib->isParameter("Time"))
00973 //      loadBasicWorksetInfo( workset, overlapped_x, overlapped_xdot, current_time );
00974       loadBasicWorksetInfo( workset, current_time );
00975     else
00976 //      loadBasicWorksetInfo( workset, overlapped_x, overlapped_xdot,
00977       loadBasicWorksetInfo( workset,
00978           paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time") );
00979 
00980     workset.params = params;
00981     workset.Vx = overlapped_Vx;
00982     workset.Vxdot = overlapped_Vxdot;
00983     workset.Vxdotdot = overlapped_Vxdotdot;
00984     workset.Vp = vp;
00985 
00986     workset.f            = overlapped_f;
00987     workset.JV           = overlapped_JV;
00988     workset.fp           = overlapped_fp;
00989     workset.j_coeff      = beta;
00990     workset.m_coeff      = alpha;
00991     workset.n_coeff      = omega;
00992 
00993     workset.num_cols_x = num_cols_x;
00994     workset.num_cols_p = num_cols_p;
00995     workset.param_offset = param_offset;
00996 
00997     workset.coord_deriv_indices = &coord_deriv_indices;
00998 
00999     for (int ws=0; ws < numWorksets; ws++) {
01000       loadWorksetBucketInfo<PHAL::AlbanyTraits::Tangent>(workset, ws);
01001       workset.ws_coord_derivs = ws_coord_derivs[ws];
01002 
01003       // FillType template argument used to specialize Sacado
01004       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Tangent>(workset);
01005       if (nfm!=Teuchos::null)
01006         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Tangent>(workset);
01007     }
01008   }
01009 
01010   vp = Teuchos::null;
01011   params = Teuchos::null;
01012 
01013   // Assemble global residual
01014   if (f != NULL)
01015     f->Export(*overlapped_f, *exporter, Add);
01016 
01017   // Assemble derivatives
01018   if (JV != NULL)
01019     JV->Export(*overlapped_JV, *exporter, Add);
01020   if (fp != NULL)
01021     fp->Export(*overlapped_fp, *exporter, Add);
01022 
01023   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
01024   if (dfm!=Teuchos::null) {
01025     PHAL::Workset workset;
01026 
01027     workset.num_cols_x = num_cols_x;
01028     workset.num_cols_p = num_cols_p;
01029     workset.param_offset = param_offset;
01030 
01031     workset.f = rcp(f,false);
01032     workset.fp = rcp(fp,false);
01033     workset.JV = rcp(JV,false);
01034     workset.j_coeff = beta;
01035     workset.n_coeff = omega;
01036     workset.x = Teuchos::rcpFromRef(x);
01037     workset.Vx = rcp(Vx,false);
01038     if (xdot != NULL) workset.transientTerms = true;
01039     if (xdotdot != NULL) workset.accelerationTerms = true;
01040 
01041     loadWorksetNodesetInfo(workset);
01042 
01043     if ( paramLib->isParameter("Time") )
01044       workset.current_time = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
01045     else
01046       workset.current_time = current_time;
01047 
01048     // FillType template argument used to specialize Sacado
01049     dfm->evaluateFields<PHAL::AlbanyTraits::Tangent>(workset);
01050   }
01051 
01052 //*out << "fp " << *fp << std::endl;
01053 
01054 }
01055 
01056 void
01057 Albany::Application::
01058 evaluateResponse(int response_index,
01059      const double current_time,
01060      const Epetra_Vector* xdot,
01061      const Epetra_Vector* xdotdot,
01062      const Epetra_Vector& x,
01063      const Teuchos::Array<ParamVec>& p,
01064      Epetra_Vector& g)
01065 {
01066   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Responses");
01067   double t = current_time;
01068   if ( paramLib->isParameter("Time") )
01069     t = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
01070 
01071   responses[response_index]->evaluateResponse(t, xdot, xdotdot, x, p, g);
01072 }
01073 
01074 void
01075 Albany::Application::
01076 evaluateResponseTangent(int response_index,
01077       const double alpha,
01078       const double beta,
01079       const double omega,
01080       const double current_time,
01081       bool sum_derivs,
01082       const Epetra_Vector* xdot,
01083       const Epetra_Vector* xdotdot,
01084       const Epetra_Vector& x,
01085       const Teuchos::Array<ParamVec>& p,
01086       ParamVec* deriv_p,
01087       const Epetra_MultiVector* Vxdot,
01088       const Epetra_MultiVector* Vxdotdot,
01089       const Epetra_MultiVector* Vx,
01090       const Epetra_MultiVector* Vp,
01091       Epetra_Vector* g,
01092       Epetra_MultiVector* gx,
01093       Epetra_MultiVector* gp)
01094 {
01095   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Response Tangent");
01096   double t = current_time;
01097   if ( paramLib->isParameter("Time") )
01098     t = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
01099 
01100   responses[response_index]->evaluateTangent(
01101     alpha, beta, omega, t, sum_derivs, xdot, xdotdot, x, p, deriv_p, Vxdot, Vxdotdot, Vx, Vp, g, gx, gp);
01102 }
01103 
01104 void
01105 Albany::Application::
01106 evaluateResponseDerivative(
01107   int response_index,
01108   const double current_time,
01109   const Epetra_Vector* xdot,
01110   const Epetra_Vector* xdotdot,
01111   const Epetra_Vector& x,
01112   const Teuchos::Array<ParamVec>& p,
01113   ParamVec* deriv_p,
01114   Epetra_Vector* g,
01115   const EpetraExt::ModelEvaluator::Derivative& dg_dx,
01116   const EpetraExt::ModelEvaluator::Derivative& dg_dxdot,
01117   const EpetraExt::ModelEvaluator::Derivative& dg_dxdotdot,
01118   const EpetraExt::ModelEvaluator::Derivative& dg_dp)
01119 {
01120   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: Response Gradient");
01121   double t = current_time;
01122   if ( paramLib->isParameter("Time") )
01123     t = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
01124 
01125   responses[response_index]->evaluateDerivative(
01126     t, xdot, xdotdot, x, p, deriv_p, g, dg_dx, dg_dxdot, dg_dxdotdot, dg_dp);
01127 }
01128 
01129 #ifdef ALBANY_SG_MP
01130 void
01131 Albany::Application::
01132 computeGlobalSGResidual(
01133   const double current_time,
01134   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
01135   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
01136   const Stokhos::EpetraVectorOrthogPoly& sg_x,
01137   const Teuchos::Array<ParamVec>& p,
01138   const Teuchos::Array<int>& sg_p_index,
01139   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
01140   Stokhos::EpetraVectorOrthogPoly& sg_f)
01141 {
01142   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: SGResidual");
01143 
01144   postRegSetup("SGResidual");
01145 
01146   //std::cout << sg_x << std::endl;
01147 
01148   // Load connectivity map and coordinates
01149   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
01150         wsElNodeEqID = disc->getWsElNodeEqID();
01151   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
01152         coords = disc->getCoords();
01153   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
01154   //      sHeight = disc->getSurfaceHeight();
01155   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
01156   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
01157 
01158   int numWorksets = wsElNodeEqID.size();
01159 
01160   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
01161   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
01162 
01163   if (sg_overlapped_x == Teuchos::null ||
01164       sg_overlapped_x->size() != sg_x.size()) {
01165     sg_overlap_map =
01166       rcp(new Epetra_LocalMap(sg_basis->size(), 0,
01167             product_comm->TimeDomainComm()));
01168     sg_overlapped_x =
01169       rcp(new Stokhos::EpetraVectorOrthogPoly(
01170       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01171     sg_overlapped_xdot =
01172   rcp(new Stokhos::EpetraVectorOrthogPoly(
01173         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01174     sg_overlapped_xdotdot =
01175   rcp(new Stokhos::EpetraVectorOrthogPoly(
01176         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01177     sg_overlapped_f =
01178       rcp(new Stokhos::EpetraVectorOrthogPoly(
01179       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01180   }
01181 
01182   for (int i=0; i<sg_x.size(); i++) {
01183 
01184     // Scatter x and xdot to the overlapped distrbution
01185     (*sg_overlapped_x)[i].Import(sg_x[i], *importer, Insert);
01186     if (sg_xdot != NULL) (*sg_overlapped_xdot)[i].Import((*sg_xdot)[i], *importer, Insert);
01187     if (sg_xdotdot != NULL) (*sg_overlapped_xdotdot)[i].Import((*sg_xdotdot)[i], *importer, Insert);
01188 
01189     // Zero out overlapped residual
01190     (*sg_overlapped_f)[i].PutScalar(0.0);
01191     sg_f[i].PutScalar(0.0);
01192 
01193   }
01194 
01195   // Set parameters
01196   for (int i=0; i<p.size(); i++)
01197     for (unsigned int j=0; j<p[i].size(); j++)
01198       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
01199 
01200   // put current_time (from Rythmos) if this is a transient problem, then compute dt
01201   //  if (sg_xdot != NULL) timeMgr.setTime(current_time);
01202 
01203 #ifdef ALBANY_CUTR
01204   if (shapeParamsHaveBeenReset) {
01205     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
01206 *out << " Calling moveMesh with params: " << std::setprecision(8);
01207 for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
01208 *out << std::endl;
01209     meshMover->moveMesh(shapeParams, morphFromInit);
01210     coords = disc->getCoords();
01211     shapeParamsHaveBeenReset = false;
01212   }
01213 #endif
01214 
01215   // Set SG parameters
01216   for (int i=0; i<sg_p_index.size(); i++) {
01217     int ii = sg_p_index[i];
01218     for (unsigned int j=0; j<p[ii].size(); j++)
01219       p[ii][j].family->setValue<PHAL::AlbanyTraits::SGResidual>(sg_p_vals[ii][j]);
01220   }
01221 
01222   // Set data in Workset struct, and perform fill via field manager
01223   {
01224     PHAL::Workset workset;
01225 
01226     workset.sg_expansion = sg_expansion;
01227     workset.sg_x         = sg_overlapped_x;
01228     workset.sg_xdot      = sg_overlapped_xdot;
01229     workset.sg_xdotdot      = sg_overlapped_xdotdot;
01230     workset.sg_f         = sg_overlapped_f;
01231 
01232     workset.current_time = current_time;
01233     //workset.delta_time = timeMgr.getDeltaTime();
01234     if (sg_xdot != NULL) workset.transientTerms = true;
01235     if (sg_xdotdot != NULL) workset.accelerationTerms = true;
01236 
01237     for (int ws=0; ws < numWorksets; ws++) {
01238       loadWorksetBucketInfo<PHAL::AlbanyTraits::SGResidual>(workset, ws);
01239 
01240       // FillType template argument used to specialize Sacado
01241       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::SGResidual>(workset);
01242       if (nfm!=Teuchos::null)
01243         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::SGResidual>(workset);
01244     }
01245   }
01246 
01247   // Assemble global residual
01248   for (int i=0; i<sg_f.size(); i++) {
01249     sg_f[i].Export((*sg_overlapped_f)[i], *exporter, Add);
01250   }
01251 
01252   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
01253   if (dfm!=Teuchos::null) {
01254     PHAL::Workset workset;
01255 
01256     workset.sg_f = Teuchos::rcpFromRef(sg_f);
01257     loadWorksetNodesetInfo(workset);
01258     workset.sg_x = Teuchos::rcpFromRef(sg_x);
01259     if (sg_xdot != NULL) workset.transientTerms = true;
01260     if (sg_xdotdot != NULL) workset.accelerationTerms = true;
01261 
01262     if ( paramLib->isParameter("Time") )
01263       workset.current_time = paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
01264     else
01265       workset.current_time = current_time;
01266 
01267     // FillType template argument used to specialize Sacado
01268     dfm->evaluateFields<PHAL::AlbanyTraits::SGResidual>(workset);
01269 
01270   }
01271 
01272   //std::cout << sg_f << std::endl;
01273 }
01274 
01275 void
01276 Albany::Application::
01277 computeGlobalSGJacobian(
01278   const double alpha,
01279   const double beta,
01280   const double omega,
01281   const double current_time,
01282   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
01283   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
01284   const Stokhos::EpetraVectorOrthogPoly& sg_x,
01285   const Teuchos::Array<ParamVec>& p,
01286   const Teuchos::Array<int>& sg_p_index,
01287   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
01288   Stokhos::EpetraVectorOrthogPoly* sg_f,
01289   Stokhos::VectorOrthogPoly<Epetra_CrsMatrix>& sg_jac)
01290 {
01291   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: SGJacobian");
01292 
01293   postRegSetup("SGJacobian");
01294 
01295   // Load connectivity map and coordinates
01296   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
01297         wsElNodeEqID = disc->getWsElNodeEqID();
01298   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
01299         coords = disc->getCoords();
01300   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
01301   //      sHeight = disc->getSurfaceHeight();
01302   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
01303   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
01304 
01305   int numWorksets = wsElNodeEqID.size();
01306 
01307   Teuchos::RCP<Epetra_CrsMatrix>& overlapped_jac = solMgr->get_overlapped_jac();
01308   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
01309   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
01310 
01311   if (sg_overlapped_x == Teuchos::null ||
01312       sg_overlapped_x->size() != sg_x.size()) {
01313     sg_overlap_map =
01314       rcp(new Epetra_LocalMap(sg_basis->size(), 0,
01315             product_comm->TimeDomainComm()));
01316     sg_overlapped_x =
01317       rcp(new Stokhos::EpetraVectorOrthogPoly(
01318       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01319     sg_overlapped_xdot =
01320   rcp(new Stokhos::EpetraVectorOrthogPoly(
01321         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01322     sg_overlapped_xdotdot =
01323   rcp(new Stokhos::EpetraVectorOrthogPoly(
01324         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01325     sg_overlapped_f =
01326       rcp(new Stokhos::EpetraVectorOrthogPoly(
01327       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01328   }
01329 
01330   for (int i=0; i<sg_x.size(); i++) {
01331 
01332     // Scatter x and xdot to the overlapped distrbution
01333     (*sg_overlapped_x)[i].Import(sg_x[i], *importer, Insert);
01334     if (sg_xdot != NULL) (*sg_overlapped_xdot)[i].Import((*sg_xdot)[i], *importer, Insert);
01335     if (sg_xdotdot != NULL) (*sg_overlapped_xdotdot)[i].Import((*sg_xdotdot)[i], *importer, Insert);
01336 
01337     // Zero out overlapped residual
01338     if (sg_f != NULL) {
01339       (*sg_overlapped_f)[i].PutScalar(0.0);
01340       (*sg_f)[i].PutScalar(0.0);
01341     }
01342 
01343   }
01344 
01345   // Create, resize and initialize overlapped Jacobians
01346   if (sg_overlapped_jac == Teuchos::null ||
01347       sg_overlapped_jac->size() != sg_jac.size()) {
01348     RCP<const Stokhos::OrthogPolyBasis<int,double> > sg_basis =
01349       sg_expansion->getBasis();
01350     RCP<Epetra_LocalMap> sg_overlap_jac_map =
01351       rcp(new Epetra_LocalMap(sg_jac.size(), 0,
01352             sg_overlap_map->Comm()));
01353     sg_overlapped_jac =
01354       rcp(new Stokhos::VectorOrthogPoly<Epetra_CrsMatrix>(
01355          sg_basis, sg_overlap_jac_map, *overlapped_jac));
01356   }
01357   for (int i=0; i<sg_overlapped_jac->size(); i++)
01358     (*sg_overlapped_jac)[i].PutScalar(0.0);
01359 
01360   // Zero out overlapped Jacobian
01361   for (int i=0; i<sg_jac.size(); i++)
01362     (*sg_overlapped_jac)[i].PutScalar(0.0);
01363 
01364   // Set parameters
01365   for (int i=0; i<p.size(); i++)
01366     for (unsigned int j=0; j<p[i].size(); j++)
01367       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
01368 
01369   // put current_time (from Rythmos) if this is a transient problem, then compute dt
01370   //  if (sg_xdot != NULL) timeMgr.setTime(current_time);
01371 
01372 #ifdef ALBANY_CUTR
01373   if (shapeParamsHaveBeenReset) {
01374     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
01375 *out << " Calling moveMesh with params: " << std::setprecision(8);
01376 for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
01377 *out << std::endl;
01378     meshMover->moveMesh(shapeParams, morphFromInit);
01379     coords = disc->getCoords();
01380     shapeParamsHaveBeenReset = false;
01381   }
01382 #endif
01383 
01384   // Set SG parameters
01385   for (int i=0; i<sg_p_index.size(); i++) {
01386     int ii = sg_p_index[i];
01387     for (unsigned int j=0; j<p[ii].size(); j++)
01388       p[ii][j].family->setValue<PHAL::AlbanyTraits::SGJacobian>(sg_p_vals[ii][j]);
01389   }
01390 
01391   RCP< Stokhos::EpetraVectorOrthogPoly > sg_overlapped_ff;
01392   if (sg_f != NULL)
01393     sg_overlapped_ff = sg_overlapped_f;
01394 
01395   // Set data in Workset struct, and perform fill via field manager
01396   {
01397     PHAL::Workset workset;
01398 
01399     workset.sg_expansion = sg_expansion;
01400     workset.sg_x         = sg_overlapped_x;
01401     workset.sg_xdot      = sg_overlapped_xdot;
01402     workset.sg_xdotdot      = sg_overlapped_xdotdot;
01403     workset.sg_f         = sg_overlapped_ff;
01404 
01405     workset.sg_Jac       = sg_overlapped_jac;
01406     loadWorksetJacobianInfo(workset, alpha, beta, omega);
01407     workset.current_time = current_time;
01408     //workset.delta_time = timeMgr.getDeltaTime();
01409     if (sg_xdot != NULL) workset.transientTerms = true;
01410     if (sg_xdotdot != NULL) workset.accelerationTerms = true;
01411 
01412     for (int ws=0; ws < numWorksets; ws++) {
01413       loadWorksetBucketInfo<PHAL::AlbanyTraits::SGJacobian>(workset, ws);
01414 
01415       // FillType template argument used to specialize Sacado
01416       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::SGJacobian>(workset);
01417       if (nfm!=Teuchos::null)
01418         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::SGJacobian>(workset);
01419     }
01420   }
01421 
01422   // Assemble global residual
01423   if (sg_f != NULL)
01424     for (int i=0; i<sg_f->size(); i++)
01425       (*sg_f)[i].Export((*sg_overlapped_f)[i], *exporter, Add);
01426 
01427   // Assemble block Jacobians
01428   RCP<Epetra_CrsMatrix> jac;
01429   for (int i=0; i<sg_jac.size(); i++) {
01430     jac = sg_jac.getCoeffPtr(i);
01431     jac->PutScalar(0.0);
01432     jac->Export((*sg_overlapped_jac)[i], *exporter, Add);
01433     jac->FillComplete(true);
01434   }
01435 
01436   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
01437   if (dfm!=Teuchos::null) {
01438     PHAL::Workset workset;
01439 
01440     workset.sg_f = rcp(sg_f,false);
01441     workset.sg_Jac = Teuchos::rcpFromRef(sg_jac);
01442     workset.j_coeff = beta;
01443     workset.n_coeff = omega;
01444     workset.sg_x = Teuchos::rcpFromRef(sg_x);
01445     if (sg_xdot != NULL) workset.transientTerms = true;
01446     if (sg_xdotdot != NULL) workset.accelerationTerms = true;
01447 
01448     loadWorksetNodesetInfo(workset);
01449 
01450     // FillType template argument used to specialize Sacado
01451     dfm->evaluateFields<PHAL::AlbanyTraits::SGJacobian>(workset);
01452   }
01453 
01454   //std::cout << sg_jac << std::endl;
01455 }
01456 
01457 void
01458 Albany::Application::
01459 computeGlobalSGTangent(
01460   const double alpha,
01461   const double beta,
01462   const double omega,
01463   const double current_time,
01464   bool sum_derivs,
01465   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
01466   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
01467   const Stokhos::EpetraVectorOrthogPoly& sg_x,
01468   const Teuchos::Array<ParamVec>& par,
01469   const Teuchos::Array<int>& sg_p_index,
01470   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
01471   ParamVec* deriv_par,
01472   const Epetra_MultiVector* Vx,
01473   const Epetra_MultiVector* Vxdot,
01474   const Epetra_MultiVector* Vxdotdot,
01475   const Epetra_MultiVector* Vp,
01476   Stokhos::EpetraVectorOrthogPoly* sg_f,
01477   Stokhos::EpetraMultiVectorOrthogPoly* sg_JVx,
01478   Stokhos::EpetraMultiVectorOrthogPoly* sg_fVp)
01479 {
01480   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: SGTangent");
01481 
01482   postRegSetup("SGTangent");
01483 
01484   // Load connectivity map and coordinates
01485   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
01486         wsElNodeEqID = disc->getWsElNodeEqID();
01487   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
01488         coords = disc->getCoords();
01489   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
01490   //      sHeight = disc->getSurfaceHeight();
01491   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
01492   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
01493 
01494   int numWorksets = wsElNodeEqID.size();
01495 
01496   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
01497   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
01498 
01499   if (sg_overlapped_x == Teuchos::null ||
01500       sg_overlapped_x->size() != sg_x.size()) {
01501     sg_overlap_map =
01502       rcp(new Epetra_LocalMap(sg_basis->size(), 0,
01503             product_comm->TimeDomainComm()));
01504     sg_overlapped_x =
01505       rcp(new Stokhos::EpetraVectorOrthogPoly(
01506       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01507     sg_overlapped_xdot =
01508   rcp(new Stokhos::EpetraVectorOrthogPoly(
01509         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01510     sg_overlapped_xdotdot =
01511   rcp(new Stokhos::EpetraVectorOrthogPoly(
01512         sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01513     sg_overlapped_f =
01514       rcp(new Stokhos::EpetraVectorOrthogPoly(
01515       sg_basis, sg_overlap_map, disc->getOverlapMap(), product_comm));
01516   }
01517 
01518   for (int i=0; i<sg_x.size(); i++) {
01519 
01520     // Scatter x and xdot to the overlapped distrbution
01521     (*sg_overlapped_x)[i].Import(sg_x[i], *importer, Insert);
01522     if (sg_xdot != NULL) (*sg_overlapped_xdot)[i].Import((*sg_xdot)[i], *importer, Insert);
01523     if (sg_xdotdot != NULL) (*sg_overlapped_xdotdot)[i].Import((*sg_xdotdot)[i], *importer, Insert);
01524 
01525     // Zero out overlapped residual
01526     if (sg_f != NULL) {
01527       (*sg_overlapped_f)[i].PutScalar(0.0);
01528       (*sg_f)[i].PutScalar(0.0);
01529     }
01530 
01531   }
01532 
01533   // Scatter Vx to the overlapped distribution
01534   RCP<Epetra_MultiVector> overlapped_Vx;
01535   if (Vx != NULL) {
01536     overlapped_Vx =
01537       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
01538          Vx->NumVectors()));
01539     overlapped_Vx->Import(*Vx, *importer, Insert);
01540   }
01541 
01542   // Scatter Vx dot to the overlapped distribution
01543   RCP<Epetra_MultiVector> overlapped_Vxdot;
01544   if (Vxdot != NULL) {
01545     overlapped_Vxdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdot->NumVectors()));
01546     overlapped_Vxdot->Import(*Vxdot, *importer, Insert);
01547   }
01548   RCP<Epetra_MultiVector> overlapped_Vxdotdot;
01549   if (Vxdotdot != NULL) {
01550     overlapped_Vxdotdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdotdot->NumVectors()));
01551     overlapped_Vxdotdot->Import(*Vxdotdot, *importer, Insert);
01552   }
01553 
01554   // Set parameters
01555   for (int i=0; i<par.size(); i++)
01556     for (unsigned int j=0; j<par[i].size(); j++)
01557       par[i][j].family->setRealValueForAllTypes(par[i][j].baseValue);
01558 
01559   // Set SG parameters
01560   for (int i=0; i<sg_p_index.size(); i++) {
01561     int ii = sg_p_index[i];
01562     for (unsigned int j=0; j<par[ii].size(); j++)
01563   par[ii][j].family->setValue<PHAL::AlbanyTraits::SGTangent>(sg_p_vals[ii][j]);
01564   }
01565 
01566   // put current_time (from Rythmos) if this is a transient problem, then compute dt
01567   //  if (sg_xdot != NULL) timeMgr.setTime(current_time);
01568 
01569   RCP<const Epetra_MultiVector > vp = rcp(Vp, false);
01570   RCP<ParamVec> params = rcp(deriv_par, false);
01571 
01572   RCP<Stokhos::EpetraVectorOrthogPoly> sg_overlapped_ff;
01573   if (sg_f != NULL)
01574     sg_overlapped_ff = sg_overlapped_f;
01575 
01576   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > sg_overlapped_JVx;
01577   if (sg_JVx != NULL) {
01578     sg_overlapped_JVx =
01579       Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
01580          sg_basis, sg_overlap_map, disc->getOverlapMap(),
01581          sg_x.productComm(),
01582          (*sg_JVx)[0].NumVectors()));
01583     sg_JVx->init(0.0);
01584   }
01585 
01586   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly > sg_overlapped_fVp;
01587   if (sg_fVp != NULL) {
01588     sg_overlapped_fVp =
01589       Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
01590          sg_basis, sg_overlap_map, disc->getOverlapMap(),
01591          sg_x.productComm(),
01592          (*sg_fVp)[0].NumVectors()));
01593     sg_fVp->init(0.0);
01594   }
01595 
01596   // Number of x & xdot tangent directions
01597   int num_cols_x = 0;
01598   if (Vx != NULL)
01599     num_cols_x = Vx->NumVectors();
01600   else if (Vxdot != NULL)
01601     num_cols_x = Vxdot->NumVectors();
01602   else if (Vxdotdot != NULL)
01603     num_cols_x = Vxdotdot->NumVectors();
01604 
01605   // Number of parameter tangent directions
01606   int num_cols_p = 0;
01607   if (params != Teuchos::null) {
01608     if (Vp != NULL)
01609       num_cols_p = Vp->NumVectors();
01610     else
01611       num_cols_p = params->size();
01612   }
01613 
01614   // Whether x and param tangent components are added or separate
01615   int param_offset = 0;
01616   if (!sum_derivs)
01617     param_offset = num_cols_x;  // offset of parameter derivs in deriv array
01618 
01619   TEUCHOS_TEST_FOR_EXCEPTION(sum_derivs &&
01620          (num_cols_x != 0) &&
01621          (num_cols_p != 0) &&
01622                      (num_cols_x != num_cols_p),
01623                      std::logic_error,
01624                      "Seed matrices Vx and Vp must have the same number " <<
01625                      " of columns when sum_derivs is true and both are "
01626                      << "non-null!" << std::endl);
01627 
01628   // Initialize
01629   if (params != Teuchos::null) {
01630     SGFadType p;
01631     int num_cols_tot = param_offset + num_cols_p;
01632     for (unsigned int i=0; i<params->size(); i++) {
01633       // Get the base value set above
01634       SGType base_val =
01635   (*params)[i].family->getValue<PHAL::AlbanyTraits::SGTangent>().val();
01636       p = SGFadType(num_cols_tot, base_val);
01637       if (Vp != NULL)
01638         for (int k=0; k<num_cols_p; k++)
01639           p.fastAccessDx(param_offset+k) = (*Vp)[k][i];
01640       else
01641         p.fastAccessDx(param_offset+i) = 1.0;
01642       (*params)[i].family->setValue<PHAL::AlbanyTraits::SGTangent>(p);
01643     }
01644   }
01645 
01646   // Set data in Workset struct, and perform fill via field manager
01647   {
01648     PHAL::Workset workset;
01649 
01650     workset.params = params;
01651     workset.sg_expansion = sg_expansion;
01652     workset.sg_x         = sg_overlapped_x;
01653     workset.sg_xdot      = sg_overlapped_xdot;
01654     workset.sg_xdotdot      = sg_overlapped_xdotdot;
01655     workset.Vx = overlapped_Vx;
01656     workset.Vxdot = overlapped_Vxdot;
01657     workset.Vxdotdot = overlapped_Vxdotdot;
01658     workset.Vp = vp;
01659 
01660     workset.sg_f         = sg_overlapped_ff;
01661     workset.sg_JV        = sg_overlapped_JVx;
01662     workset.sg_fp        = sg_overlapped_fVp;
01663     workset.j_coeff      = beta;
01664     workset.m_coeff      = alpha;
01665     workset.n_coeff      = omega;
01666 
01667     workset.num_cols_x = num_cols_x;
01668     workset.num_cols_p = num_cols_p;
01669     workset.param_offset = param_offset;
01670 
01671     workset.current_time = current_time; //timeMgr.getCurrentTime();
01672     //    workset.delta_time = timeMgr.getDeltaTime();
01673     if (sg_xdot != NULL) workset.transientTerms = true;
01674     if (sg_xdotdot != NULL) workset.accelerationTerms = true;
01675 
01676     for (int ws=0; ws < numWorksets; ws++) {
01677       loadWorksetBucketInfo<PHAL::AlbanyTraits::SGTangent>(workset, ws);
01678 
01679       // FillType template argument used to specialize Sacado
01680       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::SGTangent>(workset);
01681       if (nfm!=Teuchos::null)
01682         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::SGTangent>(workset);
01683     }
01684   }
01685 
01686   vp = Teuchos::null;
01687   params = Teuchos::null;
01688 
01689   // Assemble global residual
01690   if (sg_f != NULL)
01691     for (int i=0; i<sg_f->size(); i++)
01692       (*sg_f)[i].Export((*sg_overlapped_f)[i], *exporter, Add);
01693 
01694   // Assemble derivatives
01695   if (sg_JVx != NULL)
01696     for (int i=0; i<sg_JVx->size(); i++)
01697       (*sg_JVx)[i].Export((*sg_overlapped_JVx)[i], *exporter, Add);
01698   if (sg_fVp != NULL) {
01699     for (int i=0; i<sg_fVp->size(); i++)
01700       (*sg_fVp)[i].Export((*sg_overlapped_fVp)[i], *exporter, Add);
01701   }
01702 
01703   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
01704   if (dfm!=Teuchos::null) {
01705     PHAL::Workset workset;
01706 
01707     workset.num_cols_x = num_cols_x;
01708     workset.num_cols_p = num_cols_p;
01709     workset.param_offset = param_offset;
01710 
01711     workset.sg_f = rcp(sg_f,false);
01712     workset.sg_fp = rcp(sg_fVp,false);
01713     workset.sg_JV = rcp(sg_JVx,false);
01714     workset.j_coeff = beta;
01715     workset.n_coeff = omega;
01716     workset.sg_x = Teuchos::rcpFromRef(sg_x);
01717     workset.Vx = rcp(Vx,false);
01718     if (sg_xdot != NULL) workset.transientTerms = true;
01719     if (sg_xdotdot != NULL) workset.accelerationTerms = true;
01720 
01721     loadWorksetNodesetInfo(workset);
01722 
01723     // FillType template argument used to specialize Sacado
01724     dfm->evaluateFields<PHAL::AlbanyTraits::SGTangent>(workset);
01725   }
01726 
01727 }
01728 
01729 void
01730 Albany::Application::
01731 evaluateSGResponse(
01732   int response_index,
01733   const double curr_time,
01734   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
01735   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
01736   const Stokhos::EpetraVectorOrthogPoly& sg_x,
01737   const Teuchos::Array<ParamVec>& p,
01738   const Teuchos::Array<int>& sg_p_index,
01739   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
01740   Stokhos::EpetraVectorOrthogPoly& sg_g)
01741 {
01742   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: SGResponses");
01743 
01744   responses[response_index]->evaluateSGResponse(
01745     curr_time, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index, sg_p_vals, sg_g);
01746 }
01747 
01748 void
01749 Albany::Application::
01750 evaluateSGResponseTangent(
01751   int response_index,
01752   const double alpha,
01753   const double beta,
01754   const double omega,
01755   const double current_time,
01756   bool sum_derivs,
01757   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
01758   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
01759   const Stokhos::EpetraVectorOrthogPoly& sg_x,
01760   const Teuchos::Array<ParamVec>& p,
01761   const Teuchos::Array<int>& sg_p_index,
01762   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
01763   ParamVec* deriv_p,
01764   const Epetra_MultiVector* Vx,
01765   const Epetra_MultiVector* Vxdot,
01766   const Epetra_MultiVector* Vxdotdot,
01767   const Epetra_MultiVector* Vp,
01768   Stokhos::EpetraVectorOrthogPoly* sg_g,
01769   Stokhos::EpetraMultiVectorOrthogPoly* sg_JV,
01770   Stokhos::EpetraMultiVectorOrthogPoly* sg_gp)
01771 {
01772   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: SGResponse Tangent");
01773 
01774   responses[response_index]->evaluateSGTangent(
01775     alpha, beta, omega, current_time, sum_derivs, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index,
01776     sg_p_vals, deriv_p, Vx, Vxdot, Vxdotdot, Vp, sg_g, sg_JV, sg_gp);
01777 }
01778 
01779 void
01780 Albany::Application::
01781 evaluateSGResponseDerivative(
01782   int response_index,
01783   const double current_time,
01784   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
01785   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
01786   const Stokhos::EpetraVectorOrthogPoly& sg_x,
01787   const Teuchos::Array<ParamVec>& p,
01788   const Teuchos::Array<int>& sg_p_index,
01789   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
01790   ParamVec* deriv_p,
01791   Stokhos::EpetraVectorOrthogPoly* sg_g,
01792   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dx,
01793   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dxdot,
01794   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dxdotdot,
01795   const EpetraExt::ModelEvaluator::SGDerivative& sg_dg_dp)
01796 {
01797   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: SGResponse Gradient");
01798 
01799   responses[response_index]->evaluateSGDerivative(
01800     current_time, sg_xdot, sg_xdotdot, sg_x, p, sg_p_index, sg_p_vals, deriv_p,
01801     sg_g, sg_dg_dx, sg_dg_dxdot, sg_dg_dxdotdot, sg_dg_dp);
01802 }
01803 
01804 void
01805 Albany::Application::
01806 computeGlobalMPResidual(
01807   const double current_time,
01808   const Stokhos::ProductEpetraVector* mp_xdot,
01809   const Stokhos::ProductEpetraVector* mp_xdotdot,
01810   const Stokhos::ProductEpetraVector& mp_x,
01811   const Teuchos::Array<ParamVec>& p,
01812   const Teuchos::Array<int>& mp_p_index,
01813   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
01814   Stokhos::ProductEpetraVector& mp_f)
01815 {
01816   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: MPResidual");
01817 
01818   postRegSetup("MPResidual");
01819 
01820   // Load connectivity map and coordinates
01821   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
01822         wsElNodeEqID = disc->getWsElNodeEqID();
01823   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
01824         coords = disc->getCoords();
01825   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
01826   //      sHeight = disc->getSurfaceHeight();
01827   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
01828   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
01829 
01830   int numWorksets = wsElNodeEqID.size();
01831 
01832   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
01833   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
01834 
01835   // Create overlapped multi-point Epetra objects
01836   if (mp_overlapped_x == Teuchos::null ||
01837       mp_overlapped_x->size() != mp_x.size()) {
01838     mp_overlapped_x =
01839       rcp(new Stokhos::ProductEpetraVector(
01840       mp_x.map(), disc->getOverlapMap(), mp_x.productComm()));
01841 
01842     if (mp_xdot != NULL)
01843       mp_overlapped_xdot =
01844   rcp(new Stokhos::ProductEpetraVector(
01845         mp_xdot->map(), disc->getOverlapMap(), mp_x.productComm()));
01846 
01847     if (mp_xdotdot != NULL)
01848       mp_overlapped_xdotdot =
01849   rcp(new Stokhos::ProductEpetraVector(
01850         mp_xdotdot->map(), disc->getOverlapMap(), mp_x.productComm()));
01851 
01852   }
01853 
01854   if (mp_overlapped_f == Teuchos::null ||
01855       mp_overlapped_f->size() != mp_f.size()) {
01856     mp_overlapped_f =
01857       rcp(new Stokhos::ProductEpetraVector(
01858       mp_f.map(), disc->getOverlapMap(), mp_x.productComm()));
01859   }
01860 
01861   for (int i=0; i<mp_x.size(); i++) {
01862 
01863     // Scatter x and xdot to the overlapped distrbution
01864     (*mp_overlapped_x)[i].Import(mp_x[i], *importer, Insert);
01865     if (mp_xdot != NULL) (*mp_overlapped_xdot)[i].Import((*mp_xdot)[i], *importer, Insert);
01866     if (mp_xdotdot != NULL) (*mp_overlapped_xdotdot)[i].Import((*mp_xdotdot)[i], *importer, Insert);
01867 
01868     // Zero out overlapped residual
01869     (*mp_overlapped_f)[i].PutScalar(0.0);
01870     mp_f[i].PutScalar(0.0);
01871 
01872   }
01873 
01874   // Set parameters
01875   for (int i=0; i<p.size(); i++)
01876     for (unsigned int j=0; j<p[i].size(); j++)
01877       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
01878 
01879   // put current_time (from Rythmos) if this is a transient problem, then compute dt
01880   //  if (mp_xdot != NULL) timeMgr.setTime(current_time);
01881 
01882 #ifdef ALBANY_CUTR
01883   if (shapeParamsHaveBeenReset) {
01884     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
01885 *out << " Calling moveMesh with params: " << std::setprecision(8);
01886 for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
01887 *out << std::endl;
01888     meshMover->moveMesh(shapeParams, morphFromInit);
01889     coords = disc->getCoords();
01890     shapeParamsHaveBeenReset = false;
01891   }
01892 #endif
01893 
01894   // Set MP parameters
01895   for (int i=0; i<mp_p_index.size(); i++) {
01896     int ii = mp_p_index[i];
01897     for (unsigned int j=0; j<p[ii].size(); j++)
01898       p[ii][j].family->setValue<PHAL::AlbanyTraits::MPResidual>(mp_p_vals[ii][j]);
01899   }
01900 
01901   // Set data in Workset struct, and perform fill via field manager
01902   {
01903     PHAL::Workset workset;
01904 
01905     workset.mp_x         = mp_overlapped_x;
01906     workset.mp_xdot      = mp_overlapped_xdot;
01907     workset.mp_xdotdot      = mp_overlapped_xdotdot;
01908     workset.mp_f         = mp_overlapped_f;
01909 
01910     workset.current_time = current_time; //timeMgr.getCurrentTime();
01911     //    workset.delta_time = timeMgr.getDeltaTime();
01912     if (mp_xdot != NULL) workset.transientTerms = true;
01913     if (mp_xdotdot != NULL) workset.accelerationTerms = true;
01914 
01915     for (int ws=0; ws < numWorksets; ws++) {
01916       loadWorksetBucketInfo<PHAL::AlbanyTraits::MPResidual>(workset, ws);
01917 
01918       // FillType template argument used to specialize Sacado
01919       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::MPResidual>(workset);
01920       if (nfm!=Teuchos::null)
01921         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::MPResidual>(workset);
01922     }
01923   }
01924 
01925   // Assemble global residual
01926   for (int i=0; i<mp_f.size(); i++) {
01927     mp_f[i].Export((*mp_overlapped_f)[i], *exporter, Add);
01928   }
01929 
01930   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
01931   if (dfm!=Teuchos::null) {
01932     PHAL::Workset workset;
01933 
01934     workset.mp_f = Teuchos::rcpFromRef(mp_f);
01935     loadWorksetNodesetInfo(workset);
01936     workset.mp_x = Teuchos::rcpFromRef(mp_x);
01937     if (mp_xdot != NULL) workset.transientTerms = true;
01938     if (mp_xdotdot != NULL) workset.accelerationTerms = true;
01939 
01940     // FillType template argument used to specialize Sacado
01941     dfm->evaluateFields<PHAL::AlbanyTraits::MPResidual>(workset);
01942 
01943   }
01944 }
01945 
01946 void
01947 Albany::Application::
01948 computeGlobalMPJacobian(
01949   const double alpha,
01950   const double beta,
01951   const double omega,
01952   const double current_time,
01953   const Stokhos::ProductEpetraVector* mp_xdot,
01954   const Stokhos::ProductEpetraVector* mp_xdotdot,
01955   const Stokhos::ProductEpetraVector& mp_x,
01956   const Teuchos::Array<ParamVec>& p,
01957   const Teuchos::Array<int>& mp_p_index,
01958   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
01959   Stokhos::ProductEpetraVector* mp_f,
01960   Stokhos::ProductContainer<Epetra_CrsMatrix>& mp_jac)
01961 {
01962   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: MPJacobian");
01963 
01964   postRegSetup("MPJacobian");
01965 
01966   // Load connectivity map and coordinates
01967   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
01968         wsElNodeEqID = disc->getWsElNodeEqID();
01969   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
01970         coords = disc->getCoords();
01971   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
01972   //      sHeight = disc->getSurfaceHeight();
01973   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
01974   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
01975 
01976   int numWorksets = wsElNodeEqID.size();
01977 
01978   Teuchos::RCP<Epetra_CrsMatrix>& overlapped_jac = solMgr->get_overlapped_jac();
01979   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
01980   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
01981 
01982   // Create overlapped multi-point Epetra objects
01983   if (mp_overlapped_x == Teuchos::null ||
01984       mp_overlapped_x->size() != mp_x.size()) {
01985     mp_overlapped_x =
01986       rcp(new Stokhos::ProductEpetraVector(
01987       mp_x.map(), disc->getOverlapMap(), mp_x.productComm()));
01988 
01989     if (mp_xdot != NULL)
01990       mp_overlapped_xdot =
01991   rcp(new Stokhos::ProductEpetraVector(
01992         mp_xdot->map(), disc->getOverlapMap(), mp_x.productComm()));
01993 
01994     if (mp_xdotdot != NULL)
01995       mp_overlapped_xdotdot =
01996   rcp(new Stokhos::ProductEpetraVector(
01997         mp_xdotdot->map(), disc->getOverlapMap(), mp_x.productComm()));
01998 
01999   }
02000 
02001   if (mp_f != NULL && (mp_overlapped_f == Teuchos::null ||
02002            mp_overlapped_f->size() != mp_f->size()))
02003     mp_overlapped_f =
02004       rcp(new Stokhos::ProductEpetraVector(
02005       mp_f->map(), disc->getOverlapMap(), mp_x.productComm()));
02006 
02007   if (mp_overlapped_jac == Teuchos::null ||
02008       mp_overlapped_jac->size() != mp_jac.size())
02009     mp_overlapped_jac =
02010       rcp(new Stokhos::ProductContainer<Epetra_CrsMatrix>(
02011       mp_jac.map(), *overlapped_jac));
02012 
02013   for (int i=0; i<mp_x.size(); i++) {
02014 
02015     // Scatter x and xdot to the overlapped distrbution
02016     (*mp_overlapped_x)[i].Import(mp_x[i], *importer, Insert);
02017     if (mp_xdot != NULL) (*mp_overlapped_xdot)[i].Import((*mp_xdot)[i], *importer, Insert);
02018     if (mp_xdotdot != NULL) (*mp_overlapped_xdotdot)[i].Import((*mp_xdotdot)[i], *importer, Insert);
02019 
02020     // Zero out overlapped residual
02021     if (mp_f != NULL) {
02022       (*mp_overlapped_f)[i].PutScalar(0.0);
02023       (*mp_f)[i].PutScalar(0.0);
02024     }
02025 
02026     mp_jac[i].PutScalar(0.0);
02027     (*mp_overlapped_jac)[i].PutScalar(0.0);
02028 
02029   }
02030 
02031   // Set parameters
02032   for (int i=0; i<p.size(); i++)
02033     for (unsigned int j=0; j<p[i].size(); j++)
02034       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
02035 
02036   // put current_time (from Rythmos) if this is a transient problem, then compute dt
02037   //  if (mp_xdot != NULL) timeMgr.setTime(current_time);
02038 
02039 #ifdef ALBANY_CUTR
02040   if (shapeParamsHaveBeenReset) {
02041     TEUCHOS_FUNC_TIME_MONITOR("Albany-Cubit MeshMover");
02042 *out << " Calling moveMesh with params: " << std::setprecision(8);
02043 for (unsigned int i=0; i<shapeParams.size(); i++) *out << shapeParams[i] << "  ";
02044 *out << std::endl;
02045     meshMover->moveMesh(shapeParams, morphFromInit);
02046     coords = disc->getCoords();
02047     shapeParamsHaveBeenReset = false;
02048   }
02049 #endif
02050 
02051   // Set MP parameters
02052   for (int i=0; i<mp_p_index.size(); i++) {
02053     int ii = mp_p_index[i];
02054     for (unsigned int j=0; j<p[ii].size(); j++)
02055       p[ii][j].family->setValue<PHAL::AlbanyTraits::MPJacobian>(mp_p_vals[ii][j]);
02056   }
02057 
02058   RCP< Stokhos::ProductEpetraVector > mp_overlapped_ff;
02059   if (mp_f != NULL)
02060     mp_overlapped_ff = mp_overlapped_f;
02061 
02062   // Set data in Workset struct, and perform fill via field manager
02063   {
02064     PHAL::Workset workset;
02065 
02066     workset.mp_x         = mp_overlapped_x;
02067     workset.mp_xdot      = mp_overlapped_xdot;
02068     workset.mp_xdotdot      = mp_overlapped_xdotdot;
02069     workset.mp_f         = mp_overlapped_ff;
02070 
02071     workset.mp_Jac       = mp_overlapped_jac;
02072     loadWorksetJacobianInfo(workset, alpha, beta, omega);
02073     workset.current_time = current_time; //timeMgr.getCurrentTime();
02074     //    workset.delta_time = timeMgr.getDeltaTime();
02075     if (mp_xdot != NULL) workset.transientTerms = true;
02076     if (mp_xdotdot != NULL) workset.accelerationTerms = true;
02077 
02078     for (int ws=0; ws < numWorksets; ws++) {
02079       loadWorksetBucketInfo<PHAL::AlbanyTraits::MPJacobian>(workset, ws);
02080 
02081       // FillType template argument used to specialize Sacado
02082       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::MPJacobian>(workset);
02083       if (nfm!=Teuchos::null)
02084         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::MPJacobian>(workset);
02085     }
02086   }
02087 
02088   // Assemble global residual
02089   if (mp_f != NULL)
02090     for (int i=0; i<mp_f->size(); i++)
02091       (*mp_f)[i].Export((*mp_overlapped_f)[i], *exporter, Add);
02092 
02093   // Assemble block Jacobians
02094   RCP<Epetra_CrsMatrix> jac;
02095   for (int i=0; i<mp_jac.size(); i++) {
02096     jac = mp_jac.getCoeffPtr(i);
02097     jac->PutScalar(0.0);
02098     jac->Export((*mp_overlapped_jac)[i], *exporter, Add);
02099     jac->FillComplete(true);
02100   }
02101 
02102   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
02103   if (dfm!=Teuchos::null) {
02104     PHAL::Workset workset;
02105 
02106     workset.mp_f = rcp(mp_f,false);
02107     workset.mp_Jac = Teuchos::rcpFromRef(mp_jac);
02108     workset.j_coeff = beta;
02109     workset.n_coeff = omega;
02110     workset.mp_x = Teuchos::rcpFromRef(mp_x);;
02111     if (mp_xdot != NULL) workset.transientTerms = true;
02112     if (mp_xdotdot != NULL) workset.accelerationTerms = true;
02113 
02114     loadWorksetNodesetInfo(workset);
02115 
02116     // FillType template argument used to specialize Sacado
02117     dfm->evaluateFields<PHAL::AlbanyTraits::MPJacobian>(workset);
02118   }
02119 }
02120 
02121 void
02122 Albany::Application::
02123 computeGlobalMPTangent(
02124   const double alpha,
02125   const double beta,
02126   const double omega, 
02127   const double current_time,
02128   bool sum_derivs,
02129   const Stokhos::ProductEpetraVector* mp_xdot,
02130   const Stokhos::ProductEpetraVector* mp_xdotdot,
02131   const Stokhos::ProductEpetraVector& mp_x,
02132   const Teuchos::Array<ParamVec>& par,
02133   const Teuchos::Array<int>& mp_p_index,
02134   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
02135   ParamVec* deriv_par,
02136   const Epetra_MultiVector* Vx,
02137   const Epetra_MultiVector* Vxdot,
02138   const Epetra_MultiVector* Vxdotdot,
02139   const Epetra_MultiVector* Vp,
02140   Stokhos::ProductEpetraVector* mp_f,
02141   Stokhos::ProductEpetraMultiVector* mp_JVx,
02142   Stokhos::ProductEpetraMultiVector* mp_fVp)
02143 {
02144   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: MPTangent");
02145 
02146   postRegSetup("MPTangent");
02147 
02148   // Load connectivity map and coordinates
02149   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
02150         wsElNodeEqID = disc->getWsElNodeEqID();
02151   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
02152         coords = disc->getCoords();
02153   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
02154   //      sHeight = disc->getSurfaceHeight();
02155   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
02156   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
02157 
02158   int numWorksets = wsElNodeEqID.size();
02159 
02160   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
02161   Teuchos::RCP<Epetra_Export>& exporter = solMgr->get_exporter();
02162 
02163   // Create overlapped multi-point Epetra objects
02164   if (mp_overlapped_x == Teuchos::null ||
02165       mp_overlapped_x->size() != mp_x.size()) {
02166     mp_overlapped_x =
02167       rcp(new Stokhos::ProductEpetraVector(
02168       mp_x.map(), disc->getOverlapMap(), mp_x.productComm()));
02169 
02170     if (mp_xdot != NULL)
02171       mp_overlapped_xdot =
02172   rcp(new Stokhos::ProductEpetraVector(
02173         mp_xdot->map(), disc->getOverlapMap(), mp_x.productComm()));
02174 
02175     if (mp_xdotdot != NULL)
02176       mp_overlapped_xdotdot =
02177   rcp(new Stokhos::ProductEpetraVector(
02178         mp_xdotdot->map(), disc->getOverlapMap(), mp_x.productComm()));
02179 
02180   }
02181 
02182   if (mp_f != NULL && (mp_overlapped_f == Teuchos::null ||
02183            mp_overlapped_f->size() != mp_f->size()))
02184     mp_overlapped_f =
02185       rcp(new Stokhos::ProductEpetraVector(
02186       mp_f->map(), disc->getOverlapMap(), mp_x.productComm()));
02187 
02188   for (int i=0; i<mp_x.size(); i++) {
02189 
02190     // Scatter x and xdot to the overlapped distrbution
02191     (*mp_overlapped_x)[i].Import(mp_x[i], *importer, Insert);
02192     if (mp_xdot != NULL) (*mp_overlapped_xdot)[i].Import((*mp_xdot)[i], *importer, Insert);
02193     if (mp_xdotdot != NULL) (*mp_overlapped_xdotdot)[i].Import((*mp_xdotdot)[i], *importer, Insert);
02194 
02195     // Zero out overlapped residual
02196     if (mp_f != NULL) {
02197       (*mp_overlapped_f)[i].PutScalar(0.0);
02198       (*mp_f)[i].PutScalar(0.0);
02199     }
02200 
02201   }
02202 
02203   // Scatter Vx to the overlapped distribution
02204   RCP<Epetra_MultiVector> overlapped_Vx;
02205   if (Vx != NULL) {
02206     overlapped_Vx =
02207       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vx->NumVectors()));
02208     overlapped_Vx->Import(*Vx, *importer, Insert);
02209   }
02210 
02211   // Scatter Vx dot to the overlapped distribution
02212   RCP<Epetra_MultiVector> overlapped_Vxdot;
02213   if (Vxdot != NULL) {
02214     overlapped_Vxdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdot->NumVectors()));
02215     overlapped_Vxdot->Import(*Vxdot, *importer, Insert);
02216   }
02217   RCP<Epetra_MultiVector> overlapped_Vxdotdot;
02218   if (Vxdotdot != NULL) {
02219     overlapped_Vxdotdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdotdot->NumVectors()));
02220     overlapped_Vxdotdot->Import(*Vxdotdot, *importer, Insert);
02221   }
02222 
02223   // Set parameters
02224   for (int i=0; i<par.size(); i++)
02225     for (unsigned int j=0; j<par[i].size(); j++)
02226       par[i][j].family->setRealValueForAllTypes(par[i][j].baseValue);
02227 
02228   // Set MP parameters
02229   for (int i=0; i<mp_p_index.size(); i++) {
02230     int ii = mp_p_index[i];
02231     for (unsigned int j=0; j<par[ii].size(); j++)
02232   par[ii][j].family->setValue<PHAL::AlbanyTraits::MPTangent>(mp_p_vals[ii][j]);
02233   }
02234 
02235   // put current_time (from Rythmos) if this is a transient problem, then compute dt
02236   //  if (mp_xdot != NULL) timeMgr.setTime(current_time);
02237 
02238   RCP<const Epetra_MultiVector > vp = rcp(Vp, false);
02239   RCP<ParamVec> params = rcp(deriv_par, false);
02240 
02241   RCP< Stokhos::ProductEpetraVector > mp_overlapped_ff;
02242   if (mp_f != NULL)
02243     mp_overlapped_ff = mp_overlapped_f;
02244 
02245   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > mp_overlapped_JVx;
02246   if (mp_JVx != NULL) {
02247     mp_overlapped_JVx =
02248       Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
02249          mp_JVx->map(), disc->getOverlapMap(), mp_x.productComm(),
02250          mp_JVx->numVectors()));
02251     mp_JVx->init(0.0);
02252   }
02253 
02254   Teuchos::RCP<Stokhos::ProductEpetraMultiVector > mp_overlapped_fVp;
02255   if (mp_fVp != NULL) {
02256     mp_overlapped_fVp =
02257       Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
02258          mp_fVp->map(), disc->getOverlapMap(), mp_x.productComm(),
02259          mp_fVp->numVectors()));
02260     mp_fVp->init(0.0);
02261   }
02262 
02263   // Number of x & xdot tangent directions
02264   int num_cols_x = 0;
02265   if (Vx != NULL)
02266     num_cols_x = Vx->NumVectors();
02267   else if (Vxdot != NULL)
02268     num_cols_x = Vxdot->NumVectors();
02269   else if (Vxdotdot != NULL)
02270     num_cols_x = Vxdotdot->NumVectors();
02271 
02272   // Number of parameter tangent directions
02273   int num_cols_p = 0;
02274   if (params != Teuchos::null) {
02275     if (Vp != NULL)
02276       num_cols_p = Vp->NumVectors();
02277     else
02278       num_cols_p = params->size();
02279   }
02280 
02281   // Whether x and param tangent components are added or separate
02282   int param_offset = 0;
02283   if (!sum_derivs)
02284     param_offset = num_cols_x;  // offset of parameter derivs in deriv array
02285 
02286   TEUCHOS_TEST_FOR_EXCEPTION(sum_derivs &&
02287            (num_cols_x != 0) &&
02288            (num_cols_p != 0) &&
02289            (num_cols_x != num_cols_p),
02290            std::logic_error,
02291            "Seed matrices Vx and Vp must have the same number " <<
02292            " of columns when sum_derivs is true and both are "
02293            << "non-null!" << std::endl);
02294 
02295   // Initialize
02296   if (params != Teuchos::null) {
02297     MPFadType p;
02298     int num_cols_tot = param_offset + num_cols_p;
02299     for (unsigned int i=0; i<params->size(); i++) {
02300       // Get the base value set above
02301       MPType base_val =
02302   (*params)[i].family->getValue<PHAL::AlbanyTraits::MPTangent>().val();
02303       p = MPFadType(num_cols_tot, base_val);
02304       if (Vp != NULL)
02305         for (int k=0; k<num_cols_p; k++)
02306           p.fastAccessDx(param_offset+k) = (*Vp)[k][i];
02307       else
02308         p.fastAccessDx(param_offset+i) = 1.0;
02309       (*params)[i].family->setValue<PHAL::AlbanyTraits::MPTangent>(p);
02310     }
02311   }
02312 
02313   // Set data in Workset struct, and perform fill via field manager
02314   {
02315     PHAL::Workset workset;
02316 
02317     workset.params = params;
02318     workset.mp_x         = mp_overlapped_x;
02319     workset.mp_xdot      = mp_overlapped_xdot;
02320     workset.mp_xdotdot      = mp_overlapped_xdotdot;
02321     workset.Vx = overlapped_Vx;
02322     workset.Vxdot = overlapped_Vxdot;
02323     workset.Vxdotdot = overlapped_Vxdotdot;
02324     workset.Vp = vp;
02325 
02326     workset.mp_f         = mp_overlapped_ff;
02327     workset.mp_JV        = mp_overlapped_JVx;
02328     workset.mp_fp        = mp_overlapped_fVp;
02329     workset.j_coeff      = beta;
02330     workset.m_coeff      = alpha;
02331     workset.n_coeff      = omega;
02332 
02333     workset.num_cols_x = num_cols_x;
02334     workset.num_cols_p = num_cols_p;
02335     workset.param_offset = param_offset;
02336 
02337     workset.current_time = current_time; //timeMgr.getCurrentTime();
02338     //    workset.delta_time = timeMgr.getDeltaTime();
02339     if (mp_xdot != NULL) workset.transientTerms = true;
02340     if (mp_xdotdot != NULL) workset.accelerationTerms = true;
02341 
02342     for (int ws=0; ws < numWorksets; ws++) {
02343       loadWorksetBucketInfo<PHAL::AlbanyTraits::MPTangent>(workset, ws);
02344 
02345       // FillType template argument used to specialize Sacado
02346       fm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::MPTangent>(workset);
02347       if (nfm!=Teuchos::null)
02348         nfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::MPTangent>(workset);
02349     }
02350   }
02351 
02352   vp = Teuchos::null;
02353   params = Teuchos::null;
02354 
02355   // Assemble global residual
02356   if (mp_f != NULL)
02357     for (int i=0; i<mp_f->size(); i++)
02358       (*mp_f)[i].Export((*mp_overlapped_f)[i], *exporter, Add);
02359 
02360   // Assemble derivatives
02361   if (mp_JVx != NULL)
02362     for (int i=0; i<mp_JVx->size(); i++)
02363       (*mp_JVx)[i].Export((*mp_overlapped_JVx)[i], *exporter, Add);
02364   if (mp_fVp != NULL)
02365     for (int i=0; i<mp_fVp->size(); i++)
02366       (*mp_fVp)[i].Export((*mp_overlapped_fVp)[i], *exporter, Add);
02367 
02368   // Apply Dirichlet conditions using dfm (Dirchelt Field Manager)
02369   if (dfm!=Teuchos::null) {
02370     PHAL::Workset workset;
02371 
02372     workset.num_cols_x = num_cols_x;
02373     workset.num_cols_p = num_cols_p;
02374     workset.param_offset = param_offset;
02375 
02376     workset.mp_f = rcp(mp_f,false);
02377     workset.mp_fp = rcp(mp_fVp,false);
02378     workset.mp_JV = rcp(mp_JVx,false);
02379     workset.j_coeff = beta;
02380     workset.n_coeff = omega;
02381     workset.mp_x = Teuchos::rcpFromRef(mp_x);
02382     workset.Vx = rcp(Vx,false);
02383     if (mp_xdot != NULL) workset.transientTerms = true;
02384     if (mp_xdotdot != NULL) workset.accelerationTerms = true;
02385 
02386     loadWorksetNodesetInfo(workset);
02387 
02388     // FillType template argument used to specialize Sacado
02389     dfm->evaluateFields<PHAL::AlbanyTraits::MPTangent>(workset);
02390   }
02391 
02392 }
02393 
02394 void
02395 Albany::Application::
02396 evaluateMPResponse(
02397   int response_index,
02398   const double curr_time,
02399   const Stokhos::ProductEpetraVector* mp_xdot,
02400   const Stokhos::ProductEpetraVector* mp_xdotdot,
02401   const Stokhos::ProductEpetraVector& mp_x,
02402   const Teuchos::Array<ParamVec>& p,
02403   const Teuchos::Array<int>& mp_p_index,
02404   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
02405   Stokhos::ProductEpetraVector& mp_g)
02406 {
02407   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: MPResponses");
02408 
02409   responses[response_index]->evaluateMPResponse(
02410     curr_time, mp_xdot, mp_xdotdot, mp_x, p, mp_p_index, mp_p_vals, mp_g);
02411 }
02412 
02413 void
02414 Albany::Application::
02415 evaluateMPResponseTangent(
02416   int response_index,
02417   const double alpha,
02418   const double beta,
02419   const double omega,
02420   const double current_time,
02421   bool sum_derivs,
02422   const Stokhos::ProductEpetraVector* mp_xdot,
02423   const Stokhos::ProductEpetraVector* mp_xdotdot,
02424   const Stokhos::ProductEpetraVector& mp_x,
02425   const Teuchos::Array<ParamVec>& p,
02426   const Teuchos::Array<int>& mp_p_index,
02427   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
02428   ParamVec* deriv_p,
02429   const Epetra_MultiVector* Vx,
02430   const Epetra_MultiVector* Vxdot,
02431   const Epetra_MultiVector* Vxdotdot,
02432   const Epetra_MultiVector* Vp,
02433   Stokhos::ProductEpetraVector* mp_g,
02434   Stokhos::ProductEpetraMultiVector* mp_JV,
02435   Stokhos::ProductEpetraMultiVector* mp_gp)
02436 {
02437   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: MPResponse Tangents");
02438 
02439   responses[response_index]->evaluateMPTangent(
02440     alpha, beta, omega, current_time, sum_derivs, mp_xdot, mp_xdotdot, mp_x, p, mp_p_index,
02441     mp_p_vals, deriv_p, Vx, Vxdot, Vxdotdot, Vp, mp_g, mp_JV, mp_gp);
02442 }
02443 
02444 void
02445 Albany::Application::
02446 evaluateMPResponseDerivative(
02447   int response_index,
02448   const double current_time,
02449   const Stokhos::ProductEpetraVector* mp_xdot,
02450   const Stokhos::ProductEpetraVector* mp_xdotdot,
02451   const Stokhos::ProductEpetraVector& mp_x,
02452   const Teuchos::Array<ParamVec>& p,
02453   const Teuchos::Array<int>& mp_p_index,
02454   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
02455   ParamVec* deriv_p,
02456   Stokhos::ProductEpetraVector* mp_g,
02457   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dx,
02458   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dxdot,
02459   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dxdotdot,
02460   const EpetraExt::ModelEvaluator::MPDerivative& mp_dg_dp)
02461 {
02462   TEUCHOS_FUNC_TIME_MONITOR("> Albany Fill: MPResponse Gradient");
02463 
02464   responses[response_index]->evaluateMPDerivative(
02465     current_time, mp_xdot, mp_xdotdot, mp_x, p, mp_p_index, mp_p_vals, deriv_p,
02466     mp_g, mp_dg_dx, mp_dg_dxdot, mp_dg_dxdotdot, mp_dg_dp);
02467 }
02468 #endif //ALBANY_SG_MP
02469 
02470 void
02471 Albany::Application::
02472 evaluateStateFieldManager(const double current_time,
02473         const Epetra_Vector* xdot,
02474         const Epetra_Vector* xdotdot,
02475         const Epetra_Vector& x)
02476 {
02477   // Load connectivity map and coordinates
02478   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type&
02479         wsElNodeEqID = disc->getWsElNodeEqID();
02480   const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type&
02481         coords = disc->getCoords();
02482   //const WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type&
02483   //     sHeight = disc->getSurfaceHeight();
02484   const WorksetArray<std::string>::type& wsEBNames = disc->getWsEBNames();
02485   const WorksetArray<int>::type& wsPhysIndex = disc->getWsPhysIndex();
02486 
02487   int numWorksets = wsElNodeEqID.size();
02488 
02489   Teuchos::RCP<Epetra_Vector>& overlapped_f = solMgr->get_overlapped_f();
02490   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
02491 
02492   // visualize state field manager
02493   if (stateGraphVisDetail>0) {
02494     bool detail = false; if (stateGraphVisDetail > 1) detail=true;
02495     *out << "Phalanx writing graphviz file for graph of Residual fill (detail ="
02496    << stateGraphVisDetail << ")"<<std::endl;
02497     *out << "Process using 'dot -Tpng -O state_phalanx_graph' \n" << std::endl;
02498     for (int ps=0; ps < sfm.size(); ps++) {
02499       std::stringstream pg; pg << "state_phalanx_graph_" << ps;
02500       sfm[ps]->writeGraphvizFile<PHAL::AlbanyTraits::Residual>(pg.str(),detail,detail);
02501     }
02502     stateGraphVisDetail = -1;
02503   }
02504 
02505   // Scatter x and xdot to the overlapped distrbution
02506   solMgr->scatterX(x, xdot, xdotdot);
02507 
02508   // Set data in Workset struct
02509   PHAL::Workset workset;
02510   loadBasicWorksetInfo( workset, current_time );
02511   workset.f = overlapped_f;
02512 
02513   // Perform fill via field manager
02514   for (int ws=0; ws < numWorksets; ws++) {
02515     loadWorksetBucketInfo<PHAL::AlbanyTraits::Residual>(workset, ws);
02516     sfm[wsPhysIndex[ws]]->evaluateFields<PHAL::AlbanyTraits::Residual>(workset);
02517   }
02518 }
02519 
02520 void Albany::Application::registerShapeParameters()
02521 {
02522   int numShParams = shapeParams.size();
02523   if (shapeParamNames.size() == 0) {
02524     shapeParamNames.resize(numShParams);
02525     for (int i=0; i<numShParams; i++)
02526        shapeParamNames[i] = Albany::strint("ShapeParam",i);
02527   }
02528   Albany::DummyParameterAccessor<PHAL::AlbanyTraits::Jacobian, SPL_Traits> * dJ =
02529    new Albany::DummyParameterAccessor<PHAL::AlbanyTraits::Jacobian, SPL_Traits>();
02530   Albany::DummyParameterAccessor<PHAL::AlbanyTraits::Tangent, SPL_Traits> * dT =
02531    new Albany::DummyParameterAccessor<PHAL::AlbanyTraits::Tangent, SPL_Traits>();
02532 #ifdef ALBANY_SG_MP
02533   Albany::DummyParameterAccessor<PHAL::AlbanyTraits::SGResidual, SPL_Traits> * dSGR =
02534    new Albany::DummyParameterAccessor<PHAL::AlbanyTraits::SGResidual, SPL_Traits>();
02535   Albany::DummyParameterAccessor<PHAL::AlbanyTraits::SGJacobian, SPL_Traits> * dSGJ =
02536    new Albany::DummyParameterAccessor<PHAL::AlbanyTraits::SGJacobian, SPL_Traits>();
02537   Albany::DummyParameterAccessor<PHAL::AlbanyTraits::MPResidual, SPL_Traits> * dMPR =
02538    new Albany::DummyParameterAccessor<PHAL::AlbanyTraits::MPResidual, SPL_Traits>();
02539   Albany::DummyParameterAccessor<PHAL::AlbanyTraits::MPJacobian, SPL_Traits> * dMPJ =
02540    new Albany::DummyParameterAccessor<PHAL::AlbanyTraits::MPJacobian, SPL_Traits>();
02541 #endif //ALBANY_SG_MP
02542 
02543   // Register Parameter for Residual fill using "this->getValue" but
02544   // create dummy ones for other type that will not be used.
02545   for (int i=0; i<numShParams; i++) {
02546     *out << "Registering Shape Param " << shapeParamNames[i] << std::endl;
02547     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::Residual, SPL_Traits>
02548       (shapeParamNames[i], this, paramLib);
02549     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::Jacobian, SPL_Traits>
02550       (shapeParamNames[i], dJ, paramLib);
02551     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::Tangent, SPL_Traits>
02552       (shapeParamNames[i], dT, paramLib);
02553 #ifdef ALBANY_SG_MP
02554     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::SGResidual, SPL_Traits>
02555       (shapeParamNames[i], dSGR, paramLib);
02556     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::SGJacobian, SPL_Traits>
02557       (shapeParamNames[i], dSGJ, paramLib);
02558     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::MPResidual, SPL_Traits>
02559       (shapeParamNames[i], dMPR, paramLib);
02560     new Sacado::ParameterRegistration<PHAL::AlbanyTraits::MPJacobian, SPL_Traits>
02561       (shapeParamNames[i], dMPJ, paramLib);
02562 #endif //ALBANY_SG_MP
02563   }
02564 }
02565 
02566 PHAL::AlbanyTraits::Residual::ScalarT&
02567 Albany::Application::getValue(const std::string& name)
02568 {
02569   int index=-1;
02570   for (unsigned int i=0; i<shapeParamNames.size(); i++) {
02571     if (name == shapeParamNames[i]) index = i;
02572   }
02573   TEUCHOS_TEST_FOR_EXCEPTION(index==-1,  std::logic_error,
02574            "Error in GatherCoordinateVector::getValue, \n" <<
02575            "   Unrecognized param name: " << name << std::endl);
02576 
02577   shapeParamsHaveBeenReset = true;
02578 
02579   return shapeParams[index];
02580 }
02581 
02582 
02583 void Albany::Application::postRegSetup(std::string eval)
02584 {
02585   if (setupSet.find(eval) != setupSet.end())  return;
02586 
02587   setupSet.insert(eval);
02588 
02589   if (eval=="Residual") {
02590     for (int ps=0; ps < fm.size(); ps++)
02591       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::Residual>(eval);
02592     if (dfm!=Teuchos::null)
02593       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::Residual>(eval);
02594     if (nfm!=Teuchos::null)
02595       for (int ps=0; ps < nfm.size(); ps++)
02596         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::Residual>(eval);
02597   }
02598   else if (eval=="Jacobian") {
02599     for (int ps=0; ps < fm.size(); ps++)
02600       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::Jacobian>(eval);
02601     if (dfm!=Teuchos::null)
02602       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::Jacobian>(eval);
02603     if (nfm!=Teuchos::null)
02604       for (int ps=0; ps < nfm.size(); ps++)
02605         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::Jacobian>(eval);
02606   }
02607   else if (eval=="Tangent") {
02608     for (int ps=0; ps < fm.size(); ps++)
02609       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::Tangent>(eval);
02610     if (dfm!=Teuchos::null)
02611       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::Tangent>(eval);
02612     if (nfm!=Teuchos::null)
02613       for (int ps=0; ps < nfm.size(); ps++)
02614         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::Tangent>(eval);
02615   }
02616 #ifdef ALBANY_SG_MP
02617   else if (eval=="SGResidual") {
02618     for (int ps=0; ps < fm.size(); ps++)
02619       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::SGResidual>(eval);
02620     if (dfm!=Teuchos::null)
02621       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::SGResidual>(eval);
02622     if (nfm!=Teuchos::null)
02623       for (int ps=0; ps < nfm.size(); ps++)
02624         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::SGResidual>(eval);
02625   }
02626   else if (eval=="SGJacobian") {
02627     for (int ps=0; ps < fm.size(); ps++)
02628       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::SGJacobian>(eval);
02629     if (dfm!=Teuchos::null)
02630       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::SGJacobian>(eval);
02631     if (nfm!=Teuchos::null)
02632       for (int ps=0; ps < nfm.size(); ps++)
02633         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::SGJacobian>(eval);
02634   }
02635   else if (eval=="SGTangent") {
02636     for (int ps=0; ps < fm.size(); ps++)
02637       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::SGTangent>(eval);
02638     if (dfm!=Teuchos::null)
02639       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::SGTangent>(eval);
02640     if (nfm!=Teuchos::null)
02641       for (int ps=0; ps < nfm.size(); ps++)
02642         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::SGTangent>(eval);
02643   }
02644   else if (eval=="MPResidual") {
02645     for (int ps=0; ps < fm.size(); ps++)
02646       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::MPResidual>(eval);
02647     if (dfm!=Teuchos::null)
02648       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::MPResidual>(eval);
02649     if (nfm!=Teuchos::null)
02650       for (int ps=0; ps < nfm.size(); ps++)
02651         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::MPResidual>(eval);
02652   }
02653   else if (eval=="MPJacobian") {
02654     for (int ps=0; ps < fm.size(); ps++)
02655       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::MPJacobian>(eval);
02656     if (dfm!=Teuchos::null)
02657       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::MPJacobian>(eval);
02658     if (nfm!=Teuchos::null)
02659       for (int ps=0; ps < nfm.size(); ps++)
02660         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::MPJacobian>(eval);
02661   }
02662   else if (eval=="MPTangent") {
02663     for (int ps=0; ps < fm.size(); ps++)
02664       fm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::MPTangent>(eval);
02665     if (dfm!=Teuchos::null)
02666       dfm->postRegistrationSetupForType<PHAL::AlbanyTraits::MPTangent>(eval);
02667     if (nfm!=Teuchos::null)
02668       for (int ps=0; ps < nfm.size(); ps++)
02669         nfm[ps]->postRegistrationSetupForType<PHAL::AlbanyTraits::MPTangent>(eval);
02670   }
02671 #endif //ALBANY_SG_MP
02672   else
02673     TEUCHOS_TEST_FOR_EXCEPTION(eval!="Known Evaluation Name",  std::logic_error,
02674              "Error in setup call \n" << " Unrecognized name: " << eval << std::endl);
02675 
02676 
02677   // Write out Phalanx Graph if requested, on Proc 0, for Resid or Jacobian
02678   if (phxGraphVisDetail>0) {
02679     bool detail = false; if (phxGraphVisDetail > 1) detail=true;
02680 
02681     if (eval=="Residual") {
02682       *out << "Phalanx writing graphviz file for graph of Residual fill (detail ="
02683            << phxGraphVisDetail << ")"<<std::endl;
02684       *out << "Process using 'dot -Tpng -O phalanx_graph' \n" << std::endl;
02685       for (int ps=0; ps < fm.size(); ps++) {
02686         std::stringstream pg; pg << "phalanx_graph_" << ps;
02687         fm[ps]->writeGraphvizFile<PHAL::AlbanyTraits::Residual>(pg.str(),detail,detail);
02688       }
02689       phxGraphVisDetail = -1;
02690     }
02691     else if (eval=="Jacobian") {
02692       *out << "Phalanx writing graphviz file for graph of Jacobian fill (detail ="
02693            << phxGraphVisDetail << ")"<<std::endl;
02694       *out << "Process using 'dot -Tpng -O phalanx_graph' \n" << std::endl;
02695       for (int ps=0; ps < fm.size(); ps++) {
02696         std::stringstream pg; pg << "phalanx_graph_jac_" << ps;
02697         fm[ps]->writeGraphvizFile<PHAL::AlbanyTraits::Jacobian>(pg.str(),detail,detail);
02698       }
02699       phxGraphVisDetail = -2;
02700     }
02701   }
02702 }
02703 
02704 RCP<Epetra_Operator>
02705 Albany::Application::buildWrappedOperator(const RCP<Epetra_Operator>& Jac,
02706                                           const RCP<Epetra_Operator>& wrapInput,
02707                                           bool reorder) const
02708 {
02709   RCP<Epetra_Operator> wrappedOp = wrapInput;
02710   // if only one block just use orignal jacobian
02711   if(blockDecomp.size()==1) return (Jac);
02712 
02713   // initialize jacobian
02714   if(wrappedOp==Teuchos::null)
02715      wrappedOp = rcp(new Teko::Epetra::StridedEpetraOperator(blockDecomp,Jac));
02716   else
02717      rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->RebuildOps();
02718 
02719   // test blocked operator for correctness
02720   if(tekoParams->get("Test Blocked Operator",false)) {
02721      bool result
02722         = rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->testAgainstFullOperator(6,1e-14);
02723 
02724      *out << "Teko: Tested operator correctness:  " << (result ? "passed" : "FAILED!") << std::endl;
02725   }
02726   return wrappedOp;
02727 }
02728 
02729 void 
02730 Albany::Application::determinePiroSolver(const Teuchos::RCP<Teuchos::ParameterList>& topLevelParams){
02731 
02732   bool hasAdapt = false;
02733 
02734   const Teuchos::RCP<Teuchos::ParameterList>& problemParams =
02735     Teuchos::sublist(topLevelParams, "Problem", true);
02736 
02737   const Teuchos::RCP<Teuchos::ParameterList>& piroParams = Teuchos::sublist(topLevelParams, "Piro");
02738 
02739   if(problemParams->isSublist("Adaptation")){ // If the user has specified adaptation on input, grab the sublist
02740 
02741      hasAdapt = true;
02742 
02743   }
02744 
02745   // If not explicitly specified, determine which Piro solver to use from the problem parameters
02746   if (!piroParams->getPtr<std::string>("Solver Type")) {
02747 
02748     const std::string secondOrder = problemParams->get("Second Order", "No");
02749 
02750     TEUCHOS_TEST_FOR_EXCEPTION(
02751         secondOrder != "No" &&
02752         secondOrder != "Velocity Verlet" &&
02753         secondOrder != "Trapezoid Rule",
02754         std::logic_error,
02755         "Invalid value for Second Order: (No, Velocity Verlet, Trapezoid Rule): " <<
02756         secondOrder <<
02757         "\n");
02758 
02759     // Populate the Piro parameter list accordingly to inform the Piro solver factory
02760     std::string piroSolverToken;
02761     if (solMethod == Steady) {
02762       piroSolverToken = "NOX";
02763     } else if (solMethod == Continuation) {
02764       piroSolverToken = hasAdapt ? "LOCA Adaptive" : "LOCA";
02765     } else if (solMethod == Transient) {
02766       piroSolverToken = (secondOrder == "No") ? "Rythmos" : secondOrder;
02767     } else {
02768       // Piro cannot handle the corresponding problem
02769       piroSolverToken = "Unsupported";
02770     }
02771 
02772     piroParams->set("Solver Type", piroSolverToken);
02773 
02774   }
02775 
02776 }
02777 
02778 /*
02779 void Albany::Application::loadBasicWorksetInfo(
02780        PHAL::Workset& workset, RCP<Epetra_Vector> overlapped_x,
02781        RCP<Epetra_Vector> overlapped_xdot, double current_time)
02782 {
02783     workset.x        = overlapped_x;
02784     workset.xdot     = overlapped_xdot;
02785     workset.current_time = current_time;
02786     //workset.delta_time = delta_time;
02787     if (overlapped_xdot != Teuchos::null) workset.transientTerms = true;
02788 }
02789 */
02790 void Albany::Application::loadBasicWorksetInfo(
02791        PHAL::Workset& workset,
02792        double current_time)
02793 {
02794     workset.x        = solMgr->get_overlapped_x();
02795     workset.xdot     = solMgr->get_overlapped_xdot();
02796     workset.xdotdot     = solMgr->get_overlapped_xdotdot();
02797     workset.current_time = current_time;
02798     //workset.delta_time = delta_time;
02799     if (workset.xdot != Teuchos::null) workset.transientTerms = true;
02800     if (workset.xdotdot != Teuchos::null) workset.accelerationTerms = true;
02801 }
02802 
02803 void Albany::Application::loadWorksetJacobianInfo(PHAL::Workset& workset,
02804                                  const double& alpha, const double& beta, const double& omega)
02805 {
02806     workset.m_coeff      = alpha;
02807     workset.n_coeff      = omega;
02808     workset.j_coeff      = beta;
02809     workset.ignore_residual = ignore_residual_in_jacobian;
02810     workset.is_adjoint   = is_adjoint;
02811 }
02812 
02813 void Albany::Application::loadWorksetNodesetInfo(PHAL::Workset& workset)
02814 {
02815     workset.nodeSets = Teuchos::rcpFromRef(disc->getNodeSets());
02816     workset.nodeSetCoords = Teuchos::rcpFromRef(disc->getNodeSetCoords());
02817 
02818 }
02819 
02820 void Albany::Application::loadWorksetSidesetInfo(PHAL::Workset& workset, const int ws)
02821 {
02822 
02823     workset.sideSets = Teuchos::rcpFromRef(disc->getSideSets(ws));
02824 
02825 }
02826 
02827 void Albany::Application::setupBasicWorksetInfo(
02828   PHAL::Workset& workset,
02829   double current_time,
02830   const Epetra_Vector* xdot,
02831   const Epetra_Vector* xdotdot,
02832   const Epetra_Vector* x,
02833   const Teuchos::Array<ParamVec>& p
02834   )
02835 {
02836 
02837   Teuchos::RCP<Epetra_Vector>& overlapped_x = solMgr->get_overlapped_x();
02838   Teuchos::RCP<Epetra_Vector>& overlapped_xdot = solMgr->get_overlapped_xdot();
02839   Teuchos::RCP<Epetra_Vector>& overlapped_xdotdot = solMgr->get_overlapped_xdotdot();
02840   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
02841 
02842   // Scatter x and xdot to the overlapped distrbution
02843   solMgr->scatterX(*x, xdot, xdotdot);
02844 
02845   // Set parameters
02846   for (int i=0; i<p.size(); i++)
02847     for (unsigned int j=0; j<p[i].size(); j++)
02848       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
02849 
02850   workset.x = overlapped_x;
02851   workset.xdot = overlapped_xdot;
02852   workset.xdotdot = overlapped_xdotdot;
02853   if (!paramLib->isParameter("Time"))
02854     workset.current_time = current_time;
02855   else
02856     workset.current_time =
02857       paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
02858   if (overlapped_xdot != Teuchos::null) workset.transientTerms = true;
02859   if (overlapped_xdotdot != Teuchos::null) workset.accelerationTerms = true;
02860 
02861   // Create Teuchos::Comm from Epetra_Comm
02862   const Epetra_Comm& comm = x->Map().Comm();
02863   workset.comm = Albany::createTeuchosCommFromMpiComm(
02864                   Albany::getMpiCommFromEpetraComm(comm));
02865 
02866   workset.x_importer = importer;
02867 }
02868 
02869 #ifdef ALBANY_SG_MP
02870 void Albany::Application::setupBasicWorksetInfo(
02871   PHAL::Workset& workset,
02872   double current_time,
02873   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
02874   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
02875   const Stokhos::EpetraVectorOrthogPoly* sg_x,
02876   const Teuchos::Array<ParamVec>& p,
02877   const Teuchos::Array<int>& sg_p_index,
02878   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals)
02879 {
02880 
02881   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
02882 
02883   // Scatter x and xdot to the overlapped distrbution
02884   for (int i=0; i<sg_x->size(); i++) {
02885     (*sg_overlapped_x)[i].Import((*sg_x)[i], *importer, Insert);
02886     if (sg_xdot != NULL) (*sg_overlapped_xdot)[i].Import((*sg_xdot)[i], *importer, Insert);
02887     if (sg_xdotdot != NULL) (*sg_overlapped_xdotdot)[i].Import((*sg_xdotdot)[i], *importer, Insert);
02888   }
02889 
02890   // Set parameters
02891   for (int i=0; i<p.size(); i++)
02892     for (unsigned int j=0; j<p[i].size(); j++)
02893       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
02894 
02895   // Set SG parameters
02896   for (int i=0; i<sg_p_index.size(); i++) {
02897     int ii = sg_p_index[i];
02898     for (unsigned int j=0; j<p[ii].size(); j++) {
02899       p[ii][j].family->setValue<PHAL::AlbanyTraits::SGResidual>(sg_p_vals[ii][j]);
02900       p[ii][j].family->setValue<PHAL::AlbanyTraits::SGTangent>(sg_p_vals[ii][j]);
02901       p[ii][j].family->setValue<PHAL::AlbanyTraits::SGJacobian>(sg_p_vals[ii][j]);
02902     }
02903   }
02904 
02905   workset.sg_expansion = sg_expansion;
02906   workset.sg_x         = sg_overlapped_x;
02907   workset.sg_xdot      = sg_overlapped_xdot;
02908   workset.sg_xdotdot      = sg_overlapped_xdotdot;
02909   if (!paramLib->isParameter("Time"))
02910     workset.current_time = current_time;
02911   else
02912     workset.current_time =
02913       paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
02914   if (sg_xdot != NULL) workset.transientTerms = true;
02915   if (sg_xdotdot != NULL) workset.accelerationTerms = true;
02916 
02917   // Create Teuchos::Comm from Epetra_Comm
02918   const Epetra_Comm& comm = sg_x->coefficientMap()->Comm();
02919   workset.comm = Albany::createTeuchosCommFromMpiComm(
02920                   Albany::getMpiCommFromEpetraComm(comm));
02921 
02922   workset.x_importer = importer;
02923 }
02924 
02925 void Albany::Application::setupBasicWorksetInfo(
02926   PHAL::Workset& workset,
02927   double current_time,
02928   const Stokhos::ProductEpetraVector* mp_xdot,
02929   const Stokhos::ProductEpetraVector* mp_xdotdot,
02930   const Stokhos::ProductEpetraVector* mp_x,
02931   const Teuchos::Array<ParamVec>& p,
02932   const Teuchos::Array<int>& mp_p_index,
02933   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals)
02934 {
02935 
02936   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
02937 
02938   // Scatter x and xdot to the overlapped distrbution
02939   for (int i=0; i<mp_x->size(); i++) {
02940     (*mp_overlapped_x)[i].Import((*mp_x)[i], *importer, Insert);
02941     if (mp_xdot != NULL) (*mp_overlapped_xdot)[i].Import((*mp_xdot)[i], *importer, Insert);
02942     if (mp_xdotdot != NULL) (*mp_overlapped_xdotdot)[i].Import((*mp_xdotdot)[i], *importer, Insert);
02943   }
02944 
02945   // Set parameters
02946   for (int i=0; i<p.size(); i++)
02947     for (unsigned int j=0; j<p[i].size(); j++)
02948       p[i][j].family->setRealValueForAllTypes(p[i][j].baseValue);
02949 
02950   // Set MP parameters
02951   for (int i=0; i<mp_p_index.size(); i++) {
02952     int ii = mp_p_index[i];
02953     for (unsigned int j=0; j<p[ii].size(); j++) {
02954       p[ii][j].family->setValue<PHAL::AlbanyTraits::MPResidual>(mp_p_vals[ii][j]);
02955       p[ii][j].family->setValue<PHAL::AlbanyTraits::MPTangent>(mp_p_vals[ii][j]);
02956       p[ii][j].family->setValue<PHAL::AlbanyTraits::MPJacobian>(mp_p_vals[ii][j]);
02957     }
02958   }
02959 
02960   workset.mp_x         = mp_overlapped_x;
02961   workset.mp_xdot      = mp_overlapped_xdot;
02962   workset.mp_xdotdot      = mp_overlapped_xdotdot;
02963   if (!paramLib->isParameter("Time"))
02964     workset.current_time = current_time;
02965   else
02966     workset.current_time =
02967       paramLib->getRealValue<PHAL::AlbanyTraits::Residual>("Time");
02968   if (mp_xdot != NULL) workset.transientTerms = true;
02969   if (mp_xdotdot != NULL) workset.accelerationTerms = true;
02970 
02971   // Create Teuchos::Comm from Epetra_Comm
02972   const Epetra_Comm& comm = mp_x->coefficientMap()->Comm();
02973   workset.comm = Albany::createTeuchosCommFromMpiComm(
02974                   Albany::getMpiCommFromEpetraComm(comm));
02975 
02976   workset.x_importer = importer;
02977 }
02978 #endif //ALBANY_SG_MP
02979 
02980 void Albany::Application::setupTangentWorksetInfo(
02981   PHAL::Workset& workset,
02982   double current_time,
02983   bool sum_derivs,
02984   const Epetra_Vector* xdot,
02985   const Epetra_Vector* xdotdot,
02986   const Epetra_Vector* x,
02987   const Teuchos::Array<ParamVec>& p,
02988   ParamVec* deriv_p,
02989   const Epetra_MultiVector* Vxdot,
02990   const Epetra_MultiVector* Vxdotdot,
02991   const Epetra_MultiVector* Vx,
02992   const Epetra_MultiVector* Vp)
02993 {
02994   setupBasicWorksetInfo(workset, current_time, xdot, xdotdot, x, p);
02995 
02996   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
02997 
02998   // Scatter Vx dot the overlapped distribution
02999   RCP<Epetra_MultiVector> overlapped_Vx;
03000   if (Vx != NULL) {
03001     overlapped_Vx =
03002       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
03003             Vx->NumVectors()));
03004     overlapped_Vx->Import(*Vx, *importer, Insert);
03005   }
03006 
03007   // Scatter Vx dot the overlapped distribution
03008   RCP<Epetra_MultiVector> overlapped_Vxdot;
03009   if (Vxdot != NULL) {
03010     overlapped_Vxdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdot->NumVectors()));
03011     overlapped_Vxdot->Import(*Vxdot, *importer, Insert);
03012   }
03013   RCP<Epetra_MultiVector> overlapped_Vxdotdot;
03014   if (Vxdotdot != NULL) {
03015     overlapped_Vxdotdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdotdot->NumVectors()));
03016     overlapped_Vxdotdot->Import(*Vxdotdot, *importer, Insert);
03017   }
03018 
03019   RCP<const Epetra_MultiVector > vp = rcp(Vp, false);
03020   RCP<ParamVec> params = rcp(deriv_p, false);
03021 
03022   // Number of x & xdot tangent directions
03023   int num_cols_x = 0;
03024   if (Vx != NULL)
03025     num_cols_x = Vx->NumVectors();
03026   else if (Vxdot != NULL)
03027     num_cols_x = Vxdot->NumVectors();
03028   else if (Vxdotdot != NULL)
03029     num_cols_x = Vxdotdot->NumVectors();
03030 
03031   // Number of parameter tangent directions
03032   int num_cols_p = 0;
03033   if (params != Teuchos::null) {
03034     if (Vp != NULL)
03035       num_cols_p = Vp->NumVectors();
03036     else
03037       num_cols_p = params->size();
03038   }
03039 
03040   // Whether x and param tangent components are added or separate
03041   int param_offset = 0;
03042   if (!sum_derivs)
03043     param_offset = num_cols_x;  // offset of parameter derivs in deriv array
03044 
03045   TEUCHOS_TEST_FOR_EXCEPTION(
03046     sum_derivs &&
03047     (num_cols_x != 0) &&
03048     (num_cols_p != 0) &&
03049     (num_cols_x != num_cols_p),
03050     std::logic_error,
03051     "Seed matrices Vx and Vp must have the same number " <<
03052     " of columns when sum_derivs is true and both are "
03053     << "non-null!" << std::endl);
03054 
03055   // Initialize
03056   if (params != Teuchos::null) {
03057     TanFadType p;
03058     int num_cols_tot = param_offset + num_cols_p;
03059     for (unsigned int i=0; i<params->size(); i++) {
03060       p = TanFadType(num_cols_tot, (*params)[i].baseValue);
03061       if (Vp != NULL)
03062         for (int k=0; k<num_cols_p; k++)
03063           p.fastAccessDx(param_offset+k) = (*Vp)[k][i];
03064       else
03065         p.fastAccessDx(param_offset+i) = 1.0;
03066       (*params)[i].family->setValue<PHAL::AlbanyTraits::Tangent>(p);
03067     }
03068   }
03069 
03070   workset.params = params;
03071   workset.Vx = overlapped_Vx;
03072   workset.Vxdot = overlapped_Vxdot;
03073   workset.Vxdotdot = overlapped_Vxdotdot;
03074   workset.Vp = vp;
03075   workset.num_cols_x = num_cols_x;
03076   workset.num_cols_p = num_cols_p;
03077   workset.param_offset = param_offset;
03078 }
03079 
03080 #ifdef ALBANY_SG_MP
03081 void Albany::Application::setupTangentWorksetInfo(
03082   PHAL::Workset& workset,
03083   double current_time,
03084   bool sum_derivs,
03085   const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
03086   const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
03087   const Stokhos::EpetraVectorOrthogPoly* sg_x,
03088   const Teuchos::Array<ParamVec>& p,
03089   ParamVec* deriv_p,
03090   const Teuchos::Array<int>& sg_p_index,
03091   const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
03092   const Epetra_MultiVector* Vxdot,
03093   const Epetra_MultiVector* Vxdotdot,
03094   const Epetra_MultiVector* Vx,
03095   const Epetra_MultiVector* Vp)
03096 {
03097   setupBasicWorksetInfo(workset, current_time, sg_xdot, sg_xdotdot, sg_x, p,
03098       sg_p_index, sg_p_vals);
03099 
03100   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
03101 
03102   // Scatter Vx dot the overlapped distribution
03103   RCP<Epetra_MultiVector> overlapped_Vx;
03104   if (Vx != NULL) {
03105     overlapped_Vx =
03106       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
03107             Vx->NumVectors()));
03108     overlapped_Vx->Import(*Vx, *importer, Insert);
03109   }
03110 
03111   // Scatter Vx dot the overlapped distribution
03112   RCP<Epetra_MultiVector> overlapped_Vxdot;
03113   if (Vxdot != NULL) {
03114     overlapped_Vxdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdot->NumVectors()));
03115     overlapped_Vxdot->Import(*Vxdot, *importer, Insert);
03116   }
03117   RCP<Epetra_MultiVector> overlapped_Vxdotdot;
03118   if (Vxdotdot != NULL) {
03119     overlapped_Vxdotdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdotdot->NumVectors()));
03120     overlapped_Vxdotdot->Import(*Vxdotdot, *importer, Insert);
03121   }
03122 
03123   RCP<const Epetra_MultiVector > vp = rcp(Vp, false);
03124   RCP<ParamVec> params = rcp(deriv_p, false);
03125 
03126   // Number of x & xdot tangent directions
03127   int num_cols_x = 0;
03128   if (Vx != NULL)
03129     num_cols_x = Vx->NumVectors();
03130   else if (Vxdot != NULL)
03131     num_cols_x = Vxdot->NumVectors();
03132   else if (Vxdotdot != NULL)
03133     num_cols_x = Vxdotdot->NumVectors();
03134 
03135   // Number of parameter tangent directions
03136   int num_cols_p = 0;
03137   if (params != Teuchos::null) {
03138     if (Vp != NULL)
03139       num_cols_p = Vp->NumVectors();
03140     else
03141       num_cols_p = params->size();
03142   }
03143 
03144   // Whether x and param tangent components are added or separate
03145   int param_offset = 0;
03146   if (!sum_derivs)
03147     param_offset = num_cols_x;  // offset of parameter derivs in deriv array
03148 
03149   TEUCHOS_TEST_FOR_EXCEPTION(
03150     sum_derivs &&
03151     (num_cols_x != 0) &&
03152     (num_cols_p != 0) &&
03153     (num_cols_x != num_cols_p),
03154     std::logic_error,
03155     "Seed matrices Vx and Vp must have the same number " <<
03156     " of columns when sum_derivs is true and both are "
03157     << "non-null!" << std::endl);
03158 
03159   // Initialize
03160   if (params != Teuchos::null) {
03161     SGFadType p;
03162     int num_cols_tot = param_offset + num_cols_p;
03163     for (unsigned int i=0; i<params->size(); i++) {
03164       // Get the base value set above
03165       SGType base_val =
03166   (*params)[i].family->getValue<PHAL::AlbanyTraits::SGTangent>().val();
03167       p = SGFadType(num_cols_tot, base_val);
03168       if (Vp != NULL)
03169         for (int k=0; k<num_cols_p; k++)
03170           p.fastAccessDx(param_offset+k) = (*Vp)[k][i];
03171       else
03172         p.fastAccessDx(param_offset+i) = 1.0;
03173       (*params)[i].family->setValue<PHAL::AlbanyTraits::SGTangent>(p);
03174     }
03175   }
03176 
03177   workset.params = params;
03178   workset.Vx = overlapped_Vx;
03179   workset.Vxdot = overlapped_Vxdot;
03180   workset.Vxdotdot = overlapped_Vxdotdot;
03181   workset.Vp = vp;
03182   workset.num_cols_x = num_cols_x;
03183   workset.num_cols_p = num_cols_p;
03184   workset.param_offset = param_offset;
03185 }
03186 
03187 
03188 void Albany::Application::setupTangentWorksetInfo(
03189   PHAL::Workset& workset,
03190   double current_time,
03191   bool sum_derivs,
03192   const Stokhos::ProductEpetraVector* mp_xdot,
03193   const Stokhos::ProductEpetraVector* mp_xdotdot,
03194   const Stokhos::ProductEpetraVector* mp_x,
03195   const Teuchos::Array<ParamVec>& p,
03196   ParamVec* deriv_p,
03197   const Teuchos::Array<int>& mp_p_index,
03198   const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
03199   const Epetra_MultiVector* Vxdot,
03200   const Epetra_MultiVector* Vxdotdot,
03201   const Epetra_MultiVector* Vx,
03202   const Epetra_MultiVector* Vp)
03203 {
03204   setupBasicWorksetInfo(workset, current_time, mp_xdot, mp_xdotdot, mp_x, p,
03205       mp_p_index, mp_p_vals);
03206 
03207   Teuchos::RCP<Epetra_Import>& importer = solMgr->get_importer();
03208 
03209   // Scatter Vx dot the overlapped distribution
03210   RCP<Epetra_MultiVector> overlapped_Vx;
03211   if (Vx != NULL) {
03212     overlapped_Vx =
03213       rcp(new Epetra_MultiVector(*(disc->getOverlapMap()),
03214             Vx->NumVectors()));
03215     overlapped_Vx->Import(*Vx, *importer, Insert);
03216   }
03217 
03218   // Scatter Vx dot the overlapped distribution
03219   RCP<Epetra_MultiVector> overlapped_Vxdot;
03220   if (Vxdot != NULL) {
03221     overlapped_Vxdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdot->NumVectors()));
03222     overlapped_Vxdot->Import(*Vxdot, *importer, Insert);
03223   }
03224   RCP<Epetra_MultiVector> overlapped_Vxdotdot;
03225   if (Vxdotdot != NULL) {
03226     overlapped_Vxdotdot = rcp(new Epetra_MultiVector(*(disc->getOverlapMap()), Vxdotdot->NumVectors()));
03227     overlapped_Vxdotdot->Import(*Vxdotdot, *importer, Insert);
03228   }
03229 
03230   RCP<const Epetra_MultiVector > vp = rcp(Vp, false);
03231   RCP<ParamVec> params = rcp(deriv_p, false);
03232 
03233   // Number of x & xdot tangent directions
03234   int num_cols_x = 0;
03235   if (Vx != NULL)
03236     num_cols_x = Vx->NumVectors();
03237   else if (Vxdot != NULL)
03238     num_cols_x = Vxdot->NumVectors();
03239   else if (Vxdotdot != NULL)
03240     num_cols_x = Vxdotdot->NumVectors();
03241 
03242   // Number of parameter tangent directions
03243   int num_cols_p = 0;
03244   if (params != Teuchos::null) {
03245     if (Vp != NULL)
03246       num_cols_p = Vp->NumVectors();
03247     else
03248       num_cols_p = params->size();
03249   }
03250 
03251   // Whether x and param tangent components are added or separate
03252   int param_offset = 0;
03253   if (!sum_derivs)
03254     param_offset = num_cols_x;  // offset of parameter derivs in deriv array
03255 
03256   TEUCHOS_TEST_FOR_EXCEPTION(
03257     sum_derivs &&
03258     (num_cols_x != 0) &&
03259     (num_cols_p != 0) &&
03260     (num_cols_x != num_cols_p),
03261     std::logic_error,
03262     "Seed matrices Vx and Vp must have the same number " <<
03263     " of columns when sum_derivs is true and both are "
03264     << "non-null!" << std::endl);
03265 
03266   // Initialize
03267   if (params != Teuchos::null) {
03268     MPFadType p;
03269     int num_cols_tot = param_offset + num_cols_p;
03270     for (unsigned int i=0; i<params->size(); i++) {
03271       // Get the base value set above
03272       MPType base_val =
03273   (*params)[i].family->getValue<PHAL::AlbanyTraits::MPTangent>().val();
03274       p = MPFadType(num_cols_tot, base_val);
03275       if (Vp != NULL)
03276         for (int k=0; k<num_cols_p; k++)
03277           p.fastAccessDx(param_offset+k) = (*Vp)[k][i];
03278       else
03279         p.fastAccessDx(param_offset+i) = 1.0;
03280       (*params)[i].family->setValue<PHAL::AlbanyTraits::MPTangent>(p);
03281     }
03282   }
03283 
03284   workset.params = params;
03285   workset.Vx = overlapped_Vx;
03286   workset.Vxdot = overlapped_Vxdot;
03287   workset.Vxdotdot = overlapped_Vxdotdot;
03288   workset.Vp = vp;
03289   workset.num_cols_x = num_cols_x;
03290   workset.num_cols_p = num_cols_p;
03291   workset.param_offset = param_offset;
03292 }
03293 #endif //ALBANY_SG_MP
03294 
03295 #ifdef ALBANY_MOR
03296 Teuchos::RCP<Albany::MORFacade> Albany::Application::getMorFacade()
03297 {
03298   return morFacade;
03299 }
03300 #endif