Ifpack_ex_VectorLaplacian_FEM.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //           Galeri: Finite Element and Matrix Generation Package
00005 //                 Copyright (2006) ETHZ/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 //
00025 // Questions about Galeri? Contact Marzio Sala (marzio.sala _AT_ gmail.com)
00026 //
00027 // ************************************************************************
00028 // @HEADER
00029 
00030 #include "Galeri_ConfigDefs.h"
00031 #ifdef HAVE_MPI
00032 #include "mpi.h"
00033 #include "Epetra_MpiComm.h"
00034 #else
00035 #include "Epetra_SerialComm.h"
00036 #endif
00037 #include "Epetra_FECrsGraph.h"
00038 #include "Epetra_FEVbrMatrix.h"
00039 #include "Epetra_FEVector.h"
00040 
00041 #include "Galeri_core_Object.h"
00042 #include "Galeri_core_Workspace.h"
00043 #include "Galeri_grid_Triangle.h"
00044 #include "Galeri_grid_Loadable.h"
00045 #include "Galeri_grid_Generator.h"
00046 #include "Galeri_quadrature_Segment.h"
00047 #include "Galeri_problem_VectorLaplacian.h"
00048 #include "Galeri_viz_MEDIT.h"
00049 
00050 #include "AztecOO.h"
00051 #include "Ifpack.h"
00052 #include "Ifpack_AdditiveSchwarz.h"
00053 
00054 using namespace Galeri;
00055 
00056 class MyVectorLaplacian 
00057 {
00058   public:
00059     static inline double 
00060     getElementLHS(const double& x, const double& y, const double& z,
00061                   const int& ieq, const int& jeq,
00062                   const double& phi_i,
00063                   const double& phi_i_derx, 
00064                   const double& phi_i_dery,
00065                   const double& phi_i_derz,
00066                   const double& phi_j,
00067                   const double& phi_j_derx,
00068                   const double& phi_j_dery,
00069                   const double& phi_j_derz)
00070     {
00071       if (ieq == 0 &&  jeq == 0)
00072         return(getEpsilon(x, y, z, 0) * (phi_i_derx * phi_j_derx + 
00073                                       phi_i_dery * phi_j_dery + 
00074                                       phi_i_derz * phi_j_derz));
00075       else if (ieq == 0 && jeq == 1)
00076         return(phi_i * phi_j);
00077       else if (ieq == 1 &&  jeq == 1)
00078         return(getEpsilon(x, y, z, 1) * (phi_i_derx * phi_j_derx + 
00079                                       phi_i_dery * phi_j_dery + 
00080                                       phi_i_derz * phi_j_derz));
00081       else if (ieq == 1 && jeq == 0)
00082         return(phi_i * phi_j);
00083       else
00084         return(0.0);
00085     }
00086 
00087     static inline double
00088     getElementRHS(const double& x, const double& y, const double& z,
00089                   const int &eq, const double& phi_i)
00090 
00091     {
00092       if (eq == 0) 
00093         return((-2.0 * getEpsilon(x, y, z, 0) * exp(x) * exp(y) +
00094                 getExactSolution('f', x, y, z, 1)) * phi_i);
00095       else
00096         return(getExactSolution('f', x, y, z, 0) * phi_i);
00097     }
00098 
00099     static inline double
00100     getBoundaryValue(const double& x, const double& y, const double& z,
00101                      const int& eq)
00102     {
00103       return(getExactSolution('f', x, y, z, eq));
00104     }
00105 
00106     static inline char
00107     getBoundaryType(const int ID, const double& x, const double& y, const double& z,
00108                     const int& eq)
00109     {
00110       return('d');
00111     }
00112 
00113     static inline double 
00114     getExactSolution(const char& what, const double& x, 
00115                      const double& y, const double& z, const int eq = 0)
00116     {
00117       if (eq == 0 || eq == 1)
00118       {
00119         if (what == 'f' || what == 'x' || what =='y')
00120           return(exp(x) * exp(y));
00121         else
00122           return (0.0);
00123       }
00124       else
00125         return(0.0);
00126     }
00127 
00128     static inline double
00129     getEpsilon(const double& x, const double& y, const double& z, const int& eq)
00130     {
00131       if (eq == 0) return(1.0);
00132       else         return(1.0);
00133     }
00134 };
00135 
00136 // =========== //
00137 // main driver //
00138 // =========== //
00139 
00140 int main(int argc, char *argv[])
00141 {
00142 #ifdef HAVE_MPI
00143   MPI_Init(&argc,&argv);
00144   Epetra_MpiComm comm(MPI_COMM_WORLD);
00145 #else
00146   Epetra_SerialComm comm;
00147 #endif
00148 
00149   Galeri::core::Workspace::setNumDimensions(2);
00150 
00151   Galeri::grid::Loadable domain, boundary;
00152 
00153   int numGlobalElementsX = 4 * comm.NumProc();
00154   int numGlobalElementsY = 4;
00155 
00156   int mx = comm.NumProc();
00157   int my = 1;
00158 
00159   Galeri::grid::Generator::
00160   getSquareWithTriangles(comm, numGlobalElementsX, numGlobalElementsY,
00161                          mx, my, domain, boundary);
00162 
00163   // ============================================================ //
00164   // We are now ready to create the linear problem.               //
00165   // First, we need to define the Epetra_Map for the matrix,      //
00166   // where each grid vertex is assigned to a different            //
00167   // processor. To keep things simple, we use a linear partition. //
00168   // Then, we allocate the matrix (A), the solution vector (LHS), //
00169   // and the right-hand side (RHS).                               //
00170   // ============================================================ //
00171   
00172   int numPDEs = 2;
00173 
00174   Galeri::problem::VectorLaplacian<MyVectorLaplacian> problem(numPDEs, "Triangle");
00175 
00176   Epetra_BlockMap matrixMap(domain.getNumGlobalVertices(), numPDEs, 0, comm);
00177   Epetra_FECrsGraph matrixGraph(Copy, matrixMap, 0);
00178   problem.createGraph(domain, matrixGraph);
00179 
00180   Epetra_FEVbrMatrix A(Copy, matrixGraph);
00181   Epetra_FEVector    LHS(matrixMap);
00182   Epetra_FEVector    RHS(matrixMap);
00183 
00184   problem.integrate(domain, A, RHS);
00185 
00186   LHS.PutScalar(0.0);
00187 
00188   problem.imposeDirichletBoundaryConditions(boundary, A, RHS, LHS);
00189 
00190   // ============================================================ //
00191   // Solving the linear system is the next step, quite easy       //
00192   // because we just call AztecOO and we wait for the solution... //
00193   // ============================================================ //
00194   
00195   Epetra_LinearProblem linearProblem(&A, &LHS, &RHS);
00196   AztecOO solver(linearProblem);
00197   solver.SetAztecOption(AZ_solver, AZ_cg);
00198   solver.SetAztecOption(AZ_precond, AZ_Jacobi);
00199   solver.SetAztecOption(AZ_subdomain_solve, AZ_icc);
00200   solver.SetAztecOption(AZ_output, 16);
00201 
00202   solver.Iterate(1550, 1e-9);
00203 
00204   Epetra_MultiVector* LHSComponent = Galeri::core::Workspace::createMultiVectorComponent(LHS);
00205 
00206   for (int ieq = 0; ieq < numPDEs; ++ieq)
00207   {
00208     Galeri::core::Workspace::extractMultiVectorComponent(LHS, ieq, *LHSComponent);
00209 
00210     char fileName[80];
00211     sprintf(fileName, "sol%d", ieq);
00212     Galeri::viz::MEDIT::write(domain, fileName, *LHSComponent);
00213 
00214     problem.setEquation(ieq);
00215     problem.computeNorms(domain, *LHSComponent);
00216   }
00217   
00218 
00219   delete LHSComponent;
00220 
00221 #ifdef HAVE_MPI
00222   MPI_Finalize();
00223 #endif
00224 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:34 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3