Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_ex_VectorLaplacian_FEM.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Galeri: Finite Element and Matrix Generation Package
5 // Copyright (2002) ETHZ/Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 // ************************************************************************
40 // @HEADER
41 
42 #include "Galeri_ConfigDefs.h"
43 #ifdef HAVE_MPI
44 #include "mpi.h"
45 #include "Epetra_MpiComm.h"
46 #else
47 #include "Epetra_SerialComm.h"
48 #endif
49 #include "Epetra_FECrsGraph.h"
50 #include "Epetra_FEVbrMatrix.h"
51 #include "Epetra_FEVector.h"
52 
53 #include "Galeri_core_Object.h"
54 #include "Galeri_core_Workspace.h"
55 #include "Galeri_grid_Triangle.h"
56 #include "Galeri_grid_Loadable.h"
57 #include "Galeri_grid_Generator.h"
58 #include "Galeri_quadrature_Segment.h"
59 #include "Galeri_problem_VectorLaplacian.h"
60 #include "Galeri_viz_MEDIT.h"
61 
62 #include "AztecOO.h"
63 #include "Ifpack.h"
64 #include "Ifpack_AdditiveSchwarz.h"
65 
66 using namespace Galeri;
67 
69 {
70  public:
71  static inline double
72  getElementLHS(const double& x, const double& y, const double& z,
73  const int& ieq, const int& jeq,
74  const double& phi_i,
75  const double& phi_i_derx,
76  const double& phi_i_dery,
77  const double& phi_i_derz,
78  const double& phi_j,
79  const double& phi_j_derx,
80  const double& phi_j_dery,
81  const double& phi_j_derz)
82  {
83  if (ieq == 0 && jeq == 0)
84  return(getEpsilon(x, y, z, 0) * (phi_i_derx * phi_j_derx +
85  phi_i_dery * phi_j_dery +
86  phi_i_derz * phi_j_derz));
87  else if (ieq == 0 && jeq == 1)
88  return(phi_i * phi_j);
89  else if (ieq == 1 && jeq == 1)
90  return(getEpsilon(x, y, z, 1) * (phi_i_derx * phi_j_derx +
91  phi_i_dery * phi_j_dery +
92  phi_i_derz * phi_j_derz));
93  else if (ieq == 1 && jeq == 0)
94  return(phi_i * phi_j);
95  else
96  return(0.0);
97  }
98 
99  static inline double
100  getElementRHS(const double& x, const double& y, const double& z,
101  const int &eq, const double& phi_i)
102 
103  {
104  if (eq == 0)
105  return((-2.0 * getEpsilon(x, y, z, 0) * exp(x) * exp(y) +
106  getExactSolution('f', x, y, z, 1)) * phi_i);
107  else
108  return(getExactSolution('f', x, y, z, 0) * phi_i);
109  }
110 
111  static inline double
112  getBoundaryValue(const double& x, const double& y, const double& z,
113  const int& eq)
114  {
115  return(getExactSolution('f', x, y, z, eq));
116  }
117 
118  static inline char
119  getBoundaryType(const int ID, const double& x, const double& y, const double& z,
120  const int& eq)
121  {
122  return('d');
123  }
124 
125  static inline double
126  getExactSolution(const char& what, const double& x,
127  const double& y, const double& z, const int eq = 0)
128  {
129  if (eq == 0 || eq == 1)
130  {
131  if (what == 'f' || what == 'x' || what =='y')
132  return(exp(x) * exp(y));
133  else
134  return (0.0);
135  }
136  else
137  return(0.0);
138  }
139 
140  static inline double
141  getEpsilon(const double& x, const double& y, const double& z, const int& eq)
142  {
143  if (eq == 0) return(1.0);
144  else return(1.0);
145  }
146 };
147 
148 // =========== //
149 // main driver //
150 // =========== //
151 
152 int main(int argc, char *argv[])
153 {
154 #ifdef HAVE_MPI
155  MPI_Init(&argc,&argv);
156  Epetra_MpiComm comm(MPI_COMM_WORLD);
157 #else
158  Epetra_SerialComm comm;
159 #endif
160 
161  Galeri::core::Workspace::setNumDimensions(2);
162 
163  Galeri::grid::Loadable domain, boundary;
164 
165  int numGlobalElementsX = 4 * comm.NumProc();
166  int numGlobalElementsY = 4;
167 
168  int mx = comm.NumProc();
169  int my = 1;
170 
171  Galeri::grid::Generator::
172  getSquareWithTriangles(comm, numGlobalElementsX, numGlobalElementsY,
173  mx, my, domain, boundary);
174 
175  // ============================================================ //
176  // We are now ready to create the linear problem. //
177  // First, we need to define the Epetra_Map for the matrix, //
178  // where each grid vertex is assigned to a different //
179  // processor. To keep things simple, we use a linear partition. //
180  // Then, we allocate the matrix (A), the solution vector (LHS), //
181  // and the right-hand side (RHS). //
182  // ============================================================ //
183 
184  int numPDEs = 2;
185 
186  Galeri::problem::VectorLaplacian<MyVectorLaplacian> problem(numPDEs, "Triangle");
187 
188  Epetra_BlockMap matrixMap(domain.getNumGlobalVertices(), numPDEs, 0, comm);
189  Epetra_FECrsGraph matrixGraph(Copy, matrixMap, 0);
190  problem.createGraph(domain, matrixGraph);
191 
192  Epetra_FEVbrMatrix A(Copy, matrixGraph);
193  Epetra_FEVector LHS(matrixMap);
194  Epetra_FEVector RHS(matrixMap);
195 
196  problem.integrate(domain, A, RHS);
197 
198  LHS.PutScalar(0.0);
199 
200  problem.imposeDirichletBoundaryConditions(boundary, A, RHS, LHS);
201 
202  // ============================================================ //
203  // Solving the linear system is the next step, quite easy //
204  // because we just call AztecOO and we wait for the solution... //
205  // ============================================================ //
206 
207  Epetra_LinearProblem linearProblem(&A, &LHS, &RHS);
208  AztecOO solver(linearProblem);
209  solver.SetAztecOption(AZ_solver, AZ_cg);
210  solver.SetAztecOption(AZ_precond, AZ_Jacobi);
211  solver.SetAztecOption(AZ_subdomain_solve, AZ_icc);
212  solver.SetAztecOption(AZ_output, 16);
213 
214  solver.Iterate(1550, 1e-9);
215 
216  Epetra_MultiVector* LHSComponent = Galeri::core::Workspace::createMultiVectorComponent(LHS);
217 
218  for (int ieq = 0; ieq < numPDEs; ++ieq)
219  {
220  Galeri::core::Workspace::extractMultiVectorComponent(LHS, ieq, *LHSComponent);
221 
222  char fileName[80];
223  sprintf(fileName, "sol%d", ieq);
224  Galeri::viz::MEDIT::write(domain, fileName, *LHSComponent);
225 
226  problem.setEquation(ieq);
227  problem.computeNorms(domain, *LHSComponent);
228  }
229 
230 
231  delete LHSComponent;
232 
233 #ifdef HAVE_MPI
234  MPI_Finalize();
235 #endif
236 }
static double getElementLHS(const double &x, const double &y, const double &z, const int &ieq, const int &jeq, const double &phi_i, const double &phi_i_derx, const double &phi_i_dery, const double &phi_i_derz, const double &phi_j, const double &phi_j_derx, const double &phi_j_dery, const double &phi_j_derz)
static double getEpsilon(const double &x, const double &y, const double &z, const int &eq)
int PutScalar(double ScalarConstant)
static double getElementRHS(const double &x, const double &y, const double &z, const int &eq, const double &phi_i)
static double getExactSolution(const char &what, const double &x, const double &y, const double &z, const int eq=0)
int main(int argc, char *argv[])
int NumProc() const
static double getBoundaryValue(const double &x, const double &y, const double &z, const int &eq)
static char getBoundaryType(const int ID, const double &x, const double &y, const double &z, const int &eq)
#define RHS(a)
Definition: MatGenFD.c:60