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 "SundanceFunctional.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 "SundanceAssembler.hpp"
00042
00043
00044 using namespace Sundance;
00045 using namespace Teuchos;
00046 using namespace Playa;
00047
00048
00049 Functional::Functional(const Mesh& mesh, const Expr& integral,
00050 const Playa::VectorType<double>& vecType)
00051 : mesh_(mesh),
00052 integral_(integral),
00053 bc_(),
00054 vecType_(vecType)
00055 {
00056
00057 }
00058
00059 Functional::Functional(
00060 const Mesh& mesh,
00061 const Expr& integral,
00062 const Expr& essentialBC,
00063 const Playa::VectorType<double>& vecType)
00064 : mesh_(mesh),
00065 integral_(integral),
00066 bc_(essentialBC),
00067 vecType_(vecType)
00068 {
00069
00070 }
00071
00072
00073 LinearProblem Functional::linearVariationalProb(
00074 const Expr& var,
00075 const Expr& varEvalPts,
00076 const Expr& unk,
00077 const Expr& fixed,
00078 const Expr& fixedEvalPts) const
00079 {
00080
00081 Array<Expr> zero(unk.size());
00082 for (int i=0; i<unk.size(); i++)
00083 {
00084 Expr z = new ZeroExpr();
00085 zero[i] = z;
00086 }
00087
00088 Expr unkEvalPts = new ListExpr(zero);
00089
00090 Expr unkParams;
00091 Expr fixedParams;
00092 Expr unkParamValues;
00093 Expr fixedParamValues;
00094
00095 RCP<EquationSet> eqn
00096 = rcp(new EquationSet(integral_, bc_,
00097 tuple(var), tuple(varEvalPts),
00098 tuple(unk), tuple(unkEvalPts),
00099 unkParams, unkParamValues,
00100 tuple(fixed), tuple(fixedEvalPts)));
00101
00102 RCP<Assembler> assembler
00103 = rcp(new Assembler(mesh_, eqn, tuple(vecType_), tuple(vecType_), false));
00104
00105 return LinearProblem(assembler);
00106 }
00107
00108 NonlinearProblem Functional
00109 ::nonlinearVariationalProb(const Expr& var,
00110 const Expr& varEvalPts,
00111 const Expr& unk,
00112 const Expr& unkEvalPts,
00113 const Expr& fixed,
00114 const Expr& fixedEvalPts) const
00115 {
00116
00117 Expr unkParams;
00118 Expr fixedParams;
00119 Expr unkParamValues;
00120 Expr fixedParamValues;
00121
00122 RCP<EquationSet> eqn
00123 = rcp(new EquationSet(integral_, bc_,
00124 tuple(var), tuple(varEvalPts),
00125 tuple(unk), tuple(unkEvalPts),
00126 fixedParams, fixedParamValues,
00127 tuple(fixed), tuple(fixedEvalPts)));
00128
00129 RCP<Assembler> assembler
00130 = rcp(new Assembler(mesh_, eqn, tuple(vecType_), tuple(vecType_), false));
00131
00132 return NonlinearProblem(assembler, unkEvalPts);
00133 }
00134
00135 FunctionalEvaluator Functional::evaluator(const Expr& var,
00136 const Expr& varEvalPts,
00137 const Expr& fixed,
00138 const Expr& fixedEvalPts) const
00139 {
00140
00141 Expr unkParams;
00142 Expr fixedParams;
00143 Expr unkParamValues;
00144 Expr fixedParamValues;
00145
00146 return FunctionalEvaluator(mesh_, integral_, bc_,
00147 var,
00148 varEvalPts,
00149 fixed, fixedEvalPts,
00150 vecType_);
00151 }
00152
00153
00154 FunctionalEvaluator Functional::evaluator(const Expr& var,
00155 const Expr& varEvalPts) const
00156 {
00157 return FunctionalEvaluator(mesh_, integral_, bc_,
00158 var,
00159 varEvalPts,
00160 vecType_);
00161 }
00162
00163
00164 namespace Sundance
00165 {
00166
00167 double L2Norm(const Mesh& mesh, const CellFilter& domain,
00168 const Expr& f, const QuadratureFamily& quad,
00169 const WatchFlag& watch)
00170 {
00171 Expr I2 = Integral(domain, f*f, quad, watch);
00172
00173 return sqrt(evaluateIntegral(mesh, I2));
00174 }
00175
00176
00177 double H1Seminorm(
00178 const Mesh& mesh,
00179 const CellFilter& filter,
00180 const Expr& f,
00181 const QuadratureFamily& quad,
00182 const WatchFlag& watch)
00183 {
00184 Expr grad = gradient(mesh.spatialDim());
00185 return L2Norm(mesh, filter, grad*f, quad, watch);
00186 }
00187
00188 double H1Norm(
00189 const Mesh& mesh,
00190 const CellFilter& filter,
00191 const Expr& f,
00192 const QuadratureFamily& quad,
00193 const WatchFlag& watch)
00194 {
00195 Expr grad = gradient(mesh.spatialDim());
00196 Expr g = grad*f;
00197 Expr I2 = Integral(filter, f*f + g*g, quad, watch);
00198
00199 return sqrt(evaluateIntegral(mesh, I2));
00200 }
00201 }