Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
65
66using 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
152int main(int argc, char *argv[])
153{
154#ifdef HAVE_MPI
155 MPI_Init(&argc,&argv);
156 Epetra_MpiComm comm(MPI_COMM_WORLD);
157#else
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}
Copy
int main(int argc, char *argv[])
#define RHS(a)
Definition MatGenFD.c:60
int PutScalar(double ScalarConstant)
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)
static double getExactSolution(const char &what, const double &x, const double &y, const double &z, const int eq=0)
static double getElementRHS(const double &x, const double &y, const double &z, const int &eq, const double &phi_i)
static char getBoundaryType(const int ID, const double &x, const double &y, const double &z, const int &eq)
static double getBoundaryValue(const double &x, const double &y, const double &z, const int &eq)