42 #include "Galeri_ConfigDefs.h"
45 #include "Epetra_MpiComm.h"
47 #include "Epetra_SerialComm.h"
49 #include "Epetra_FECrsGraph.h"
50 #include "Epetra_FEVbrMatrix.h"
51 #include "Epetra_FEVector.h"
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"
66 using namespace Galeri;
73 const int& ieq,
const int& jeq,
75 const double& phi_i_derx,
76 const double& phi_i_dery,
77 const double& phi_i_derz,
79 const double& phi_j_derx,
80 const double& phi_j_dery,
81 const double& phi_j_derz)
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);
101 const int &eq,
const double& phi_i)
105 return((-2.0 * getEpsilon(x, y, z, 0) * exp(x) * exp(y) +
106 getExactSolution(
'f', x, y, z, 1)) * phi_i);
108 return(getExactSolution(
'f', x, y, z, 0) * phi_i);
115 return(getExactSolution(
'f', x, y, z, eq));
127 const double& y,
const double& z,
const int eq = 0)
129 if (eq == 0 || eq == 1)
131 if (what ==
'f' || what ==
'x' || what ==
'y')
132 return(exp(x) * exp(y));
141 getEpsilon(
const double& x,
const double& y,
const double& z,
const int& eq)
143 if (eq == 0)
return(1.0);
152 int main(
int argc,
char *argv[])
155 MPI_Init(&argc,&argv);
161 Galeri::core::Workspace::setNumDimensions(2);
163 Galeri::grid::Loadable domain, boundary;
165 int numGlobalElementsX = 4 * comm.
NumProc();
166 int numGlobalElementsY = 4;
171 Galeri::grid::Generator::
172 getSquareWithTriangles(comm, numGlobalElementsX, numGlobalElementsY,
173 mx, my, domain, boundary);
186 Galeri::problem::VectorLaplacian<MyVectorLaplacian> problem(numPDEs,
"Triangle");
188 Epetra_BlockMap matrixMap(domain.getNumGlobalVertices(), numPDEs, 0, comm);
190 problem.createGraph(domain, matrixGraph);
196 problem.integrate(domain, A, RHS);
200 problem.imposeDirichletBoundaryConditions(boundary, A, RHS, LHS);
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);
214 solver.Iterate(1550, 1e-9);
216 Epetra_MultiVector* LHSComponent = Galeri::core::Workspace::createMultiVectorComponent(LHS);
218 for (
int ieq = 0; ieq < numPDEs; ++ieq)
220 Galeri::core::Workspace::extractMultiVectorComponent(LHS, ieq, *LHSComponent);
223 sprintf(fileName,
"sol%d", ieq);
224 Galeri::viz::MEDIT::write(domain, fileName, *LHSComponent);
226 problem.setEquation(ieq);
227 problem.computeNorms(domain, *LHSComponent);
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[])
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)