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 "SundanceLinearEigenproblem.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceTestFunction.hpp"
00035 #include "SundanceUnknownFunction.hpp"
00036 #include "SundanceEssentialBC.hpp"
00037 #include "SundanceIntegral.hpp"
00038 #include "SundanceListExpr.hpp"
00039 #include "SundanceZeroExpr.hpp"
00040 #include "SundanceEquationSet.hpp"
00041 #include "SundanceQuadratureFamily.hpp"
00042 #include "SundanceAssembler.hpp"
00043 #include "SundanceMaximalCellFilter.hpp"
00044 #include "SundanceGaussianQuadrature.hpp"
00045
00046 #include "PlayaLinearCombinationDecl.hpp"
00047 #include "PlayaLinearCombinationImpl.hpp"
00048 #include "PlayaLinearOperatorDecl.hpp"
00049 #include "PlayaVectorDecl.hpp"
00050 #include "PlayaSimpleDiagonalOpDecl.hpp"
00051
00052 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00053 #include "PlayaVectorImpl.hpp"
00054 #include "PlayaLinearCombinationImpl.hpp"
00055 #include "PlayaLinearOperatorImpl.hpp"
00056 #include "PlayaSimpleDiagonalOpImpl.hpp"
00057 #endif
00058
00059 using namespace Sundance;
00060 using namespace Sundance;
00061 using namespace Sundance;
00062 using namespace Sundance;
00063 using namespace Teuchos;
00064 using namespace Playa;
00065 using namespace PlayaExprTemplates;
00066
00067 static Time& normalizationTimer()
00068 {
00069 static RCP<Time> rtn
00070 = TimeMonitor::getNewTimer("Eigenfunction normalization");
00071 return *rtn;
00072 }
00073
00074 static Time& makeEigensystemTimer()
00075 {
00076 static RCP<Time> rtn
00077 = TimeMonitor::getNewTimer("Building eigensystem stiffness matrix");
00078 return *rtn;
00079 }
00080
00081 LinearEigenproblem::LinearEigenproblem(
00082 const Mesh& mesh, const Expr& eqn,
00083 const Expr& v, const Expr& u,
00084 const VectorType<double>& vecType)
00085 : lumpMass_(false),
00086 kProb_(),
00087 mProb_(),
00088 M_(),
00089 MUnlumped_(),
00090 discSpace_()
00091 {
00092 Expr empty;
00093
00094 kProb_ = LinearProblem(mesh, eqn, empty, v, u, vecType);
00095 discSpace_ = *(kProb_.solnSpace()[0]);
00096 }
00097
00098 LinearEigenproblem::LinearEigenproblem(
00099 const Mesh& mesh, const Expr& eqn,
00100 const Expr& v, const Expr& u,
00101 const VectorType<double>& vecType,
00102 bool lumpedMass)
00103 : lumpMass_(lumpedMass),
00104 kProb_(),
00105 mProb_(),
00106 M_(),
00107 MUnlumped_(),
00108 discSpace_()
00109 {
00110 Expr empty;
00111
00112 kProb_ = LinearProblem(mesh, eqn, empty, v, u, vecType);
00113 mProb_ = makeMassProb(mesh, empty, v, u, vecType);
00114 discSpace_ = *(kProb_.solnSpace()[0]);
00115 MUnlumped_ = mProb_.getOperator();
00116 if (lumpMass_)
00117 {
00118 M_ = lumpedOperator(MUnlumped_);
00119 }
00120 else
00121 {
00122 M_ = MUnlumped_;
00123 }
00124 }
00125
00126
00127 LinearEigenproblem::LinearEigenproblem(
00128 const Mesh& mesh, const Expr& eqn,
00129 const Expr& massExpr,
00130 const Expr& v, const Expr& u,
00131 const VectorType<double>& vecType,
00132 bool lumpedMass)
00133 : lumpMass_(lumpedMass),
00134 kProb_(),
00135 mProb_(),
00136 M_(),
00137 MUnlumped_(),
00138 discSpace_()
00139 {
00140 Expr bc;
00141 kProb_ = LinearProblem(mesh, eqn, bc, v, u, vecType);
00142 mProb_ = makeMassProb(mesh, massExpr, v, u, vecType);
00143 discSpace_ = *(kProb_.solnSpace()[0]);
00144
00145 MUnlumped_ = mProb_.getOperator();
00146 if (lumpMass_)
00147 {
00148 M_ = lumpedOperator(MUnlumped_);
00149 }
00150 else
00151 {
00152 M_ = MUnlumped_;
00153 }
00154
00155 }
00156
00157 LinearProblem LinearEigenproblem::makeMassProb(
00158 const Mesh& mesh,
00159 const Expr& massExpr,
00160 const Expr& v, const Expr& u,
00161 const VectorType<double>& vecType) const
00162 {
00163 Expr eqn;
00164
00165 CellFilter interior = new MaximalCellFilter();
00166 QuadratureFamily quad = new GaussianQuadrature( 4 );
00167 if (massExpr.ptr().get()==0)
00168 {
00169 eqn = Integral(interior, v*u, quad);
00170 }
00171 else
00172 {
00173 eqn = Integral(interior, massExpr, quad);
00174 }
00175 Expr bc;
00176 LinearProblem rtn(mesh, eqn, bc, v, u, vecType);
00177 return rtn;
00178 }
00179
00180
00181 Array<Expr> LinearEigenproblem::makeEigenfunctions(
00182 Array<Vector<double> >& ev) const
00183 {
00184 TimeMonitor timer(normalizationTimer());
00185
00186 Array<Expr> x(ev.size());
00187 CellFilter interior = new MaximalCellFilter();
00188 QuadratureFamily q = new GaussianQuadrature(2);
00189 for (int i=0; i<ev.size(); i++)
00190 {
00191 x[i] = new DiscreteFunction(discSpace_, ev[i], "ev[" + Teuchos::toString(i)+"]");
00192 double N = 1.0;
00193 if (MUnlumped_.ptr().get())
00194 {
00195 N = ev[i] * (MUnlumped_ * ev[i]);
00196 }
00197 else
00198 {
00199 N = evaluateIntegral(discSpace_.mesh(),
00200 Integral(interior, x[i]*x[i], q));
00201 }
00202 ev[i].scale(1.0/sqrt(N));
00203 }
00204
00205 return x;
00206 }
00207
00208
00209 LinearOperator<double>
00210 LinearEigenproblem::lumpedOperator(const LinearOperator<double>& M) const
00211 {
00212 Vector<double> ones = M.domain().createMember();
00213 ones.setToConstant(1.0);
00214 Vector<double> m = M * ones;
00215 LinearOperator<double> rtn = diagonalOperator(m);
00216
00217 return rtn;
00218 }
00219
00220
00221 Eigensolution LinearEigenproblem::solve(const Eigensolver<double>& solver) const
00222 {
00223 Array<std::complex<double> > ew;
00224 Array<Vector<double> > ev;
00225
00226 LinearOperator<double> K;
00227 {
00228 TimeMonitor timer(makeEigensystemTimer());
00229 K = kProb_.getOperator();
00230 }
00231
00232 solver.solve(K, M_, ev, ew);
00233
00234 Array<Expr> eigenfuncs = makeEigenfunctions(ev);
00235
00236 return Eigensolution(eigenfuncs, ew);
00237 }
00238