00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceL2Projector.hpp"
00032 #include "PlayaAztecSolver.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "PlayaTabs.hpp"
00035
00036 #include "SundanceDerivative.hpp"
00037 #include "SundanceTestFunction.hpp"
00038 #include "SundanceUnknownFunction.hpp"
00039 #include "SundanceEssentialBC.hpp"
00040 #include "SundanceIntegral.hpp"
00041 #include "SundanceGaussianQuadrature.hpp"
00042
00043 #include "SundanceCellFilter.hpp"
00044 #include "SundanceMaximalCellFilter.hpp"
00045
00046 using namespace Sundance;
00047 using namespace Teuchos;
00048 using namespace Playa;
00049
00050
00051 L2Projector::L2Projector(const DiscreteSpace& space,
00052 const Expr& expr,
00053 const LinearSolver<double>& solver)
00054 : prob_(), solver_()
00055 {
00056 CoordinateSystem cs = new CartesianCoordinateSystem();
00057 init(space, cs, expr, solver, new GaussianQuadrature(4));
00058 }
00059
00060
00061 L2Projector::L2Projector(const DiscreteSpace& space,
00062 const CoordinateSystem& cs,
00063 const Expr& expr,
00064 const LinearSolver<double>& solver)
00065 : prob_(), solver_()
00066 {
00067 init(space, cs, expr, solver, new GaussianQuadrature(4));
00068 }
00069
00070
00071 L2Projector::L2Projector(const DiscreteSpace& space,
00072 const Expr& expr,
00073 const LinearSolver<double>& solver,
00074 const QuadratureFamily& quad)
00075 : prob_(), solver_()
00076 {
00077 CoordinateSystem cs = new CartesianCoordinateSystem();
00078 init(space, cs, expr, solver, quad);
00079 }
00080
00081 L2Projector::L2Projector(const DiscreteSpace& space,
00082 const Expr& expr)
00083 : prob_(), solver_()
00084 {
00085 CoordinateSystem cs = new CartesianCoordinateSystem();
00086
00087
00088 std::map<int,int> azOptions;
00089 std::map<int,double> azParams;
00090
00091 azOptions[AZ_solver] = AZ_cg;
00092 azOptions[AZ_precond] = AZ_dom_decomp;
00093 azOptions[AZ_subdomain_solve] = AZ_icc;
00094 azOptions[AZ_graph_fill] = 1;
00095 azOptions[AZ_max_iter] = 1000;
00096 azOptions[AZ_output] = AZ_none;
00097 azParams[AZ_tol] = 1.0e-13;
00098
00099 LinearSolver<double> solver = new AztecSolver(azOptions,azParams);
00100
00101 init(space, cs, expr, solver, new GaussianQuadrature(4));
00102 }
00103
00104 L2Projector::L2Projector(const DiscreteSpace& space,
00105 const Expr& expr,
00106 const QuadratureFamily& quad)
00107 : prob_(), solver_()
00108 {
00109 CoordinateSystem cs = new CartesianCoordinateSystem();
00110
00111
00112 std::map<int,int> azOptions;
00113 std::map<int,double> azParams;
00114
00115 azOptions[AZ_solver] = AZ_cg;
00116 azOptions[AZ_precond] = AZ_dom_decomp;
00117 azOptions[AZ_subdomain_solve] = AZ_icc;
00118 azOptions[AZ_graph_fill] = 1;
00119 azOptions[AZ_max_iter] = 1000;
00120 azOptions[AZ_output] = AZ_none;
00121 azParams[AZ_tol] = 1.0e-13;
00122
00123 LinearSolver<double> solver = new AztecSolver(azOptions,azParams);
00124
00125 init(space, cs, expr, solver, quad);
00126 }
00127
00128
00129 L2Projector::L2Projector(const DiscreteSpace& space,
00130 const CoordinateSystem& cs,
00131 const Expr& expr)
00132 : prob_(), solver_()
00133 {
00134
00135 std::map<int,int> azOptions;
00136 std::map<int,double> azParams;
00137
00138 azOptions[AZ_solver] = AZ_cg;
00139 azOptions[AZ_precond] = AZ_dom_decomp;
00140 azOptions[AZ_subdomain_solve] = AZ_icc;
00141 azOptions[AZ_graph_fill] = 1;
00142 azOptions[AZ_max_iter] = 1000;
00143 azOptions[AZ_output] = AZ_none;
00144 azParams[AZ_tol] = 1.0e-13;
00145
00146 LinearSolver<double> solver = new AztecSolver(azOptions,azParams);
00147
00148 init(space, cs, expr, solver, new GaussianQuadrature(4));
00149 }
00150
00151
00152
00153
00154 void L2Projector::init(const DiscreteSpace& space,
00155 const CoordinateSystem& coordSys,
00156 const Expr& expr,
00157 const LinearSolver<double>& solver,
00158 const QuadratureFamily& quad)
00159 {
00160 TEUCHOS_TEST_FOR_EXCEPTION(space.basis().size() != expr.size(),
00161 std::runtime_error,
00162 "mismatched vector structure between basis and expr");
00163
00164 TEUCHOS_TEST_FOR_EXCEPTION(space.basis().size() == 0,
00165 std::runtime_error,
00166 "Empty basis?");
00167
00168 Expr v = new TestFunction(space.basis()[0], "dummy_v[0]");
00169 Expr u = new UnknownFunction(space.basis()[0], "dummy_u[0]");
00170
00171 for (int i=1; i<space.basis().size(); i++)
00172 {
00173 v.append(new TestFunction(space.basis()[i], "dummy_v["
00174 + Teuchos::toString(i)+"]"));
00175 u.append(new UnknownFunction(space.basis()[i], "dummy_u["
00176 + Teuchos::toString(i)+"]"));
00177 }
00178
00179 CellFilter interior = new MaximalCellFilter();
00180
00181 Expr eqn = 0.0;
00182 Expr J = coordSys.jacobian();
00183
00184 for (int i=0; i<space.basis().size(); i++)
00185 {
00186 eqn = eqn + Integral(space.cellFilters(i),
00187 J*v[i]*(u[i]-expr[i]),
00188 quad);
00189 }
00190 Expr bc;
00191
00192 prob_ = LinearProblem(space.mesh(), eqn, bc, v, u, space.vecType());
00193 solver_ = solver;
00194 }
00195
00196
00197