epetra/1Dfem/ExampleApplication.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                     Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #include "ExampleApplication.hpp"
00030 #ifdef HAVE_MPI
00031 #  include "Epetra_MpiComm.h"
00032 #  include "mpi.h"
00033 #else
00034 #  include "Epetra_SerialComm.h"
00035 #endif // HAVE_MPI
00036 #include "Epetra_CrsMatrix.h"
00037 
00038 #ifdef EXAMPLEAPPLICATION_DEBUG
00039 #include <iostream>
00040 #endif
00041 
00042 #include "Teuchos_dyn_cast.hpp"
00043 
00044 ExampleApplication1Dfem::ExampleApplication1Dfem(Teuchos::RefCountPtr<Epetra_Comm> &epetra_comm_ptr, Teuchos::ParameterList &params)
00045 {
00046   initialize(epetra_comm_ptr, params);
00047 }
00048 
00049 void ExampleApplication1Dfem::initialize(Teuchos::RefCountPtr<Epetra_Comm> &epetra_comm_ptr, Teuchos::ParameterList &params)
00050 {
00051   epetra_comm_ptr_ = epetra_comm_ptr;
00052   numElements_ = params.get<int>( "NumElements" );
00053   
00054   // This object is derived from NOX::Epetra::Interface
00055   problemInterfacePtr_ = Teuchos::rcp(new TransientInterface(numElements_, *epetra_comm_ptr_, -20.0, 20.0));
00056 
00057   // Set the PDE nonlinear coefficient for this problem
00058   problemInterfacePtr_->setPDEfactor(1.0);
00059 
00060   // This is needed to extract the Epetra_Map for the solution
00061   epetra_map_ptr_ = Teuchos::rcp( new Epetra_Map( problemInterfacePtr_->getMap() ) );
00062 //  Epetra_Vector& soln = problemInterfacePtr_->getSolution();
00063 //  const Epetra_BlockMap& solnBlockMap = soln.Map();
00064 //  std::cout << "typeid(solnBlockMap).name() = " << typeid(solnBlockMap).name() << std::endl;
00065 //  const Epetra_Map& solnMap = Teuchos::dyn_cast<const Epetra_Map>(solnBlockMap);
00066 //  epetra_map_ptr_ = Teuchos::rcp( new Epetra_Map( solnMap ) );
00067 
00068   // This is needed to extract the Epetra_CrsGraph for the Jacobian
00069 //  Epetra_CrsMatrix &jacobian = problemInterfacePtr_->getJacobian();
00070 //  W_graph_ = Teuchos::rcp(new Epetra_CrsGraph( jacobian.Graph() ) );
00071   W_graph_ = Teuchos::rcp(new Epetra_CrsGraph( problemInterfacePtr_->getGraph() ) );
00072 
00073 }
00074 
00075 // Overridden from EpetraExt::ModelEvaluator
00076 
00077 Teuchos::RefCountPtr<const Epetra_Map>
00078 ExampleApplication1Dfem::get_x_map() const
00079 {
00080   return epetra_map_ptr_;
00081 }
00082 
00083 Teuchos::RefCountPtr<const Epetra_Map>
00084 ExampleApplication1Dfem::get_f_map() const
00085 {
00086   return epetra_map_ptr_;
00087 }
00088 
00089 Teuchos::RefCountPtr<const Epetra_Vector>
00090 ExampleApplication1Dfem::get_x_init() const
00091 {
00092   Epetra_Vector& soln = problemInterfacePtr_->getSolution();
00093   Teuchos::RefCountPtr<Epetra_Vector> x_init = Teuchos::rcp(new Epetra_Vector(soln));
00094   return x_init;
00095 }
00096 
00097 Teuchos::RefCountPtr<const Epetra_Vector>
00098 ExampleApplication1Dfem::get_x_dot_init() const
00099 {
00100   Epetra_Vector& soln = problemInterfacePtr_->getSolution();
00101   Teuchos::RefCountPtr<Epetra_Vector> x_dot_init = Teuchos::rcp(new Epetra_Vector(soln));
00102   x_dot_init->PutScalar(0.0);
00103   return x_dot_init;
00104 }
00105 
00106 Teuchos::RefCountPtr<Epetra_Operator>
00107 ExampleApplication1Dfem::create_W() const
00108 {
00109   Teuchos::RefCountPtr<Epetra_Operator> W = Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00110   return W;
00111 }
00112 
00113 EpetraExt::ModelEvaluator::InArgs
00114 ExampleApplication1Dfem::createInArgs() const
00115 {
00116   InArgsSetup inArgs;
00117   inArgs.setSupports(IN_ARG_x,true);
00118   inArgs.setSupports(IN_ARG_x_dot,true);
00119   inArgs.setSupports(IN_ARG_alpha,true);
00120   inArgs.setSupports(IN_ARG_beta,true);
00121   return inArgs;
00122 }
00123 
00124 EpetraExt::ModelEvaluator::OutArgs
00125 ExampleApplication1Dfem::createOutArgs() const
00126 {
00127   OutArgsSetup outArgs;
00128   outArgs.setSupports(OUT_ARG_f,true);
00129   outArgs.setSupports(OUT_ARG_W,true);
00130   return outArgs;
00131 }
00132 
00133 void ExampleApplication1Dfem::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
00134 {
00135   Teuchos::RefCountPtr<const Epetra_Vector> x = inArgs.get_x();
00136   Teuchos::RefCountPtr<const Epetra_Vector> xdot = inArgs.get_x_dot();
00137 #ifdef EXAMPLEAPPLICATION_DEBUG
00138   std::cout << "ExampleApplication1Dfem::evalModel ---------------------------{" << std::endl;
00139   std::cout << "x = " << std::endl;
00140   x->Print(std::cout);
00141   std::cout << "xdot = " << std::endl;
00142   xdot->Print(std::cout);
00143 #endif // EXAMPLEAPPLICATION_DEBUG
00144   Teuchos::RefCountPtr<Epetra_Vector> f;
00145   if( (f = outArgs.get_f()).get() ) 
00146   {
00147     NOX::Epetra::Interface::Required::FillType flag = NOX::Epetra::Interface::Required::Residual;
00148     problemInterfacePtr_->evaluate(flag,&*x,&*xdot,0.0,0.0,&*f,NULL);
00149 #ifdef EXAMPLEAPPLICATION_DEBUG
00150     std::cout << "f = " << std::endl;
00151     f->Print(std::cout);
00152 #endif // EXAMPLEAPPLICATION_DEBUG
00153   }
00154   Teuchos::RefCountPtr<Epetra_Operator> W;
00155   if( (W = outArgs.get_W()).get() ) 
00156   {
00157     const double alpha = inArgs.get_alpha();
00158     const double beta = inArgs.get_beta();
00159     NOX::Epetra::Interface::Required::FillType flag = NOX::Epetra::Interface::Required::Jac;
00160 //    Epetra_CrsMatrix& jacobian = problemInterfacePtr_->getJacobian();
00161     Epetra_CrsMatrix& jac = Teuchos::dyn_cast<Epetra_CrsMatrix>(*W);
00162     problemInterfacePtr_->evaluate(flag,&*x,&*xdot,alpha,beta,NULL,&jac);
00163 #ifdef EXAMPLEAPPLICATION_DEBUG
00164     std::cout << "jac = " << std::endl;
00165     jac.Print(std::cout);
00166 #endif // EXAMPLEAPPLICATION_DEBUG
00167   }
00168 #ifdef EXAMPLEAPPLICATION_DEBUG
00169   std::cout << "ExampleApplication1Dfem::evalModel ---------------------------}" << std::endl;
00170 #endif // EXAMPLEAPPLICATION_DEBUG
00171 }
00172 
00173 Teuchos::RefCountPtr<Epetra_Vector> ExampleApplication1Dfem::get_exact_solution( double t ) const
00174 {
00175   Epetra_Vector& x_exact = problemInterfacePtr_->getExactSoln(t);
00176   Teuchos::RefCountPtr<Epetra_Vector> x_exact_ptr = Teuchos::rcp(&x_exact,false);
00177   return(x_exact_ptr);
00178 }

Generated on Thu Sep 18 12:30:05 2008 for Rythmos - Transient Integration for Differential Equations by doxygen 1.3.9.1