00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <SolverExpression.hpp>
00021
00022 #include <Solver.hpp>
00023 #include <PDEProblem.hpp>
00024 #include <PDESystem.hpp>
00025 #include <SolverDriver.hpp>
00026 #include <BoundaryConditionSet.hpp>
00027
00028 #include <FunctionExpressionValue.hpp>
00029
00030 #include <PDESystemExpression.hpp>
00031 #include <VariationalProblemExpression.hpp>
00032
00033 #include <ScalarDiscretizationTypeBase.hpp>
00034 #include <ScalarDiscretizationTypeSpectral.hpp>
00035
00036 #include <DiscretizationType.hpp>
00037
00038 #include <ScalarFunctionBase.hpp>
00039
00040 #include <FEMSolution.hpp>
00041 #include <LegendreSolution.hpp>
00042
00043 #include <VariationalProblem.hpp>
00044
00045 #include <SceneExpression.hpp>
00046 #include <Scene.hpp>
00047
00048 #include <Structured3DMesh.hpp>
00049 #include <Information.hpp>
00050
00051 #include <DegreeOfFreedomSetBuilder.hpp>
00052
00053 #include <ParameterCenter.hpp>
00054
00055 #include <SpectralFunction.hpp>
00056 #include <FEMFunction.hpp>
00057 #include <FEMFunctionBuilder.hpp>
00058
00059 #include <DomainExpression.hpp>
00060 #include <DomainExpressionSet.hpp>
00061 #include <DomainExpressionVariable.hpp>
00062 #include <DomainExpressionAnalytic.hpp>
00063
00064 #include <VariableRepository.hpp>
00065
00066
00067 std::ostream& SolverExpression::
00068 put(std::ostream& os) const
00069 {
00070 os << __FILE__ << ':' << __LINE__ << ": NOT IMPLEMENTED\n";
00071 return os;
00072 }
00073
00074
00075 SolverExpression::
00076 SolverExpression(ReferenceCounting<UnknownListExpression> unknownList,
00077 ReferenceCounting<MeshExpression> mesh,
00078 ReferenceCounting<SolverOptionsExpression> solverOptions,
00079 ReferenceCounting<ProblemExpression> problemExpression,
00080 ReferenceCounting<DomainExpression> domainExpression)
00081 : Expression(Expression::solver),
00082 __unknownList(unknownList),
00083 __mesh(mesh),
00084 __solverOptions(solverOptions),
00085 __problemExpression(problemExpression),
00086 __domain(domainExpression)
00087 {
00088 ;
00089 }
00090
00091 SolverExpression::
00092 ~SolverExpression()
00093 {
00094 ;
00095 }
00096
00097 void SolverExpression::
00098 __setScene(ReferenceCounting<DomainExpression> domainExp) const
00099 {
00100 switch(domainExp->domainType()) {
00101 case DomainExpression::set: {
00102 DomainExpressionSet& D
00103 = dynamic_cast<DomainExpressionSet&>(*domainExp);
00104
00105 const SceneExpression& SE
00106 = dynamic_cast<const SceneExpression&>(*(D.scene()));
00107 Information::instance().setScene(SE.scene());
00108 break;
00109 }
00110 case DomainExpression::analytic: {
00111 DomainExpressionAnalytic& D
00112 = dynamic_cast<DomainExpressionAnalytic&>(*domainExp);
00113
00114 Information::instance().setScene(D.domain()->scene());
00115 break;
00116 }
00117 case DomainExpression::variable: {
00118 DomainExpressionVariable& D
00119 = dynamic_cast<DomainExpressionVariable&>(*domainExp);
00120 __setScene(D.domainExpression());
00121 break;
00122 }
00123 case DomainExpression::undefined: {
00124 throw ErrorHandler(__FILE__,__LINE__,
00125 "undefined domain",
00126 ErrorHandler::normal);
00127 }
00128 default: {
00129 throw ErrorHandler(__FILE__,__LINE__,
00130 "unexpected domain type",
00131 ErrorHandler::unexpected);
00132 }
00133 }
00134 }
00135
00136 void SolverExpression::
00137 execute()
00138 {
00140 Information::instance().setUnknownList(__unknownList);
00141
00143 ParameterCenter::instance().reset();
00144
00145 __mesh->execute();
00146
00148 Information::instance().setMesh(__mesh->mesh());
00149
00150 __unknownList->execute();
00151
00152 if (__domain != 0) {
00153 __domain->execute();
00154 __setScene(__domain);
00155 }
00156
00157 __solverOptions->execute();
00158 __problemExpression->execute();
00159
00160 if (__domain != 0) {
00161 switch (__problemExpression->problemType()) {
00162 case ProblemExpression::pdeSystem: {
00163 (*dynamic_cast<PDESystemExpression&>(*__problemExpression).pdeSystem()).setDomain(__domain->domain());
00164 break;
00165 }
00166 case ProblemExpression::variationalProblem: {
00167 (*dynamic_cast<VariationalProblemExpression&>(*__problemExpression).variationalProblem()).setDomain(__domain->domain());
00168 break;
00169 }
00170 default: {
00171 throw ErrorHandler(__FILE__,__LINE__,
00172 "unexpected problem type",
00173 ErrorHandler::unexpected);
00174 }
00175 }
00176 } else {
00177 if (__problemExpression->hasPOVBoundary()) {
00178 throw ErrorHandler(__FILE__,__LINE__,
00179 "cannot use POVRay references in standard FEM discretization.\n"
00180 "You probably forgot to specify the domain in the 'solve' bloc.\n"
00181 "Refere to the Fictitious Domain part of the documentation",
00182 ErrorHandler::normal);
00183 }
00184 if (__problemExpression->hasPredefinedBoundary()) {
00185 throw ErrorHandler(__FILE__,__LINE__,
00186 "cannot external surface mesh in standard FEM discretization.\n"
00187 "You probably forgot to specify the domain in the 'solve' bloc.\n"
00188 "Refere to the Fictitious Domain part of the documentation",
00189 ErrorHandler::normal);
00190 }
00191 }
00192
00193 DiscretizationType discretization;
00194
00195 for (UnknownListExpression::iterator i= __unknownList->begin();
00196 i != __unknownList->end(); ++i) {
00197 UnknownExpression& unknown = **i;
00198 unknown.execute();
00199
00200 ReferenceCounting<ScalarDiscretizationTypeBase> unknownDiscretization = 0;
00201
00202 ScalarDiscretizationTypeBase::Type unknownDiscretizationType
00203 = ScalarDiscretizationTypeBase::getDefault(unknown.discretizationType());
00204
00205 switch(unknownDiscretizationType) {
00206 case ScalarDiscretizationTypeBase::functionLike: {
00207 ReferenceCounting<FunctionVariable> unknownVariable
00208 = VariableRepository::instance().findVariable<FunctionVariable>(unknown.name());
00209
00210 ConstReferenceCounting<ScalarFunctionBase> f = unknownVariable->expression()->function();
00211 switch (f->type()) {
00212 case ScalarFunctionBase::femfunction: {
00213 unknownDiscretization
00214 = new ScalarDiscretizationTypeFEM(dynamic_cast<const FEMFunctionBase&>(*f).discretizationType());
00215 break;
00216 }
00217 case ScalarFunctionBase::spectral: {
00218 const SpectralFunction& spectralFunction
00219 = dynamic_cast<const SpectralFunction&>(*f);
00220 unknownDiscretization
00221 = new ScalarDiscretizationTypeSpectral(spectralFunction.mesh()->degrees(),
00222 spectralFunction.mesh()->shape().a(),
00223 spectralFunction.mesh()->shape().b());
00224 break;
00225 }
00226 default: {
00227 unknownDiscretization
00228 = new ScalarDiscretizationTypeFEM(ScalarDiscretizationTypeBase::lagrangianFEM1);
00229 }
00230 }
00231 unknown.setDiscretizationType(unknownDiscretization->type());
00232 break;
00233 }
00234 case ScalarDiscretizationTypeBase::spectralLegendre: {
00235 ReferenceCounting<FunctionVariable> unknownVariable
00236 = VariableRepository::instance().findVariable<FunctionVariable>(unknown.name());
00237
00238 ConstReferenceCounting<ScalarFunctionBase> f = unknownVariable->expression()->function();
00239 switch (f->type()) {
00240 case ScalarFunctionBase::spectral: {
00241 const SpectralFunction& spectralFunction
00242 = dynamic_cast<const SpectralFunction&>(*f);
00243 unknownDiscretization
00244 = new ScalarDiscretizationTypeSpectral(spectralFunction.mesh()->degrees(),
00245 spectralFunction.mesh()->shape().a(),
00246 spectralFunction.mesh()->shape().b());
00247 break;
00248 }
00249 default: {
00250 throw ErrorHandler(__FILE__,__LINE__,
00251 "unknown is not a spectral function",
00252 ErrorHandler::unexpected);
00253 }
00254 }
00255 unknown.setDiscretizationType(unknownDiscretization->type());
00256 break;
00257 }
00258 case ScalarDiscretizationTypeBase::lagrangianFEM0:
00259 case ScalarDiscretizationTypeBase::lagrangianFEM1:
00260 case ScalarDiscretizationTypeBase::lagrangianFEM2: {
00261 unknownDiscretization
00262 = new ScalarDiscretizationTypeFEM(unknownDiscretizationType);
00263 break;
00264 }
00265 case ScalarDiscretizationTypeBase::undefined: {
00266
00267 unknownDiscretization
00268 = new ScalarDiscretizationTypeFEM(ScalarDiscretizationTypeBase::lagrangianFEM1);
00269 }
00270 }
00271
00272 ASSERT(unknownDiscretization != 0);
00273
00274 discretization.add(unknownDiscretization);
00275 }
00276
00277 ASSERT(discretization.number()>0);
00278
00279 switch (discretization.getFamily()) {
00280 case DiscretizationType::fem: {
00281 this->__solveFEM(discretization);
00282 break;
00283 }
00284 case DiscretizationType::spectral: {
00285 this->__solveLegendre(discretization);
00286 break;
00287 }
00288 default: {
00289 throw ErrorHandler(__FILE__,__LINE__,
00290 "Discretization '("+stringify(discretization)+")' is not possible",
00291 ErrorHandler::normal);
00292 }
00293 }
00294
00296 Information::instance().unsetUnknownList();
00297
00299 Information::instance().unsetMesh();
00300 }
00301
00302 void SolverExpression::
00303 __solveFEM(const DiscretizationType& discretization)
00304 {
00305 ReferenceCounting<DegreeOfFreedomSetBuilder> dofBuilder;
00306 if (__domain != 0) {
00307 dofBuilder = new DegreeOfFreedomSetBuilder(discretization,
00308 *__mesh->mesh(),
00309 *__domain->domain());
00310 } else {
00311 dofBuilder = new DegreeOfFreedomSetBuilder(discretization,
00312 *__mesh->mesh());
00313 }
00314
00315 const DegreeOfFreedomSet& degreeOfFreedomSet
00316 = dofBuilder->degreeOfFreedomSet();
00317
00318 FEMSolution u(degreeOfFreedomSet);
00319
00320 typedef
00321 std::vector<ReferenceCounting<FEMFunctionBase> >
00322 UnknownList;
00323 UnknownList UL;
00324 ConstReferenceCounting<Mesh> mesh = __mesh->mesh();
00325
00326 for (UnknownListExpression::iterator i= __unknownList->begin();
00327 i != __unknownList->end(); ++i) {
00328 UnknownExpression& unknown = **i;
00329
00330 const size_t unknownNumber = __unknownList->number(unknown.name());
00331
00332 ReferenceCounting<FunctionVariable> unknownVariable
00333 = VariableRepository::instance().findVariable<FunctionVariable>(unknown.name());
00334
00335 ConstReferenceCounting<ScalarFunctionBase> f = unknownVariable->expression()->function();
00336
00337 FEMFunctionBuilder femFunctionBuilder;
00338 femFunctionBuilder.build(dynamic_cast<const ScalarDiscretizationTypeFEM&>(discretization[unknownNumber]), mesh, *f);
00339 UL.push_back(femFunctionBuilder.getBuiltFEMFunction());
00340 }
00341
00342 u.setUserFunction(UL);
00343
00344 ConstReferenceCounting<Problem> P;
00345 switch (__problemExpression->problemType()) {
00346 case ProblemExpression::pdeSystem: {
00347 P = (PDESystem*)(dynamic_cast<PDESystemExpression&>(*__problemExpression).pdeSystem());
00348 break;
00349 }
00350 case ProblemExpression::variationalProblem: {
00351 P = (VariationalProblem*)(dynamic_cast<VariationalProblemExpression&>(*__problemExpression).variationalProblem());
00352 break;
00353 }
00354 default: {
00355 throw ErrorHandler(__FILE__,__LINE__,
00356 "unexpected problem type",
00357 ErrorHandler::unexpected);
00358 }
00359 }
00360
00361 ffout(2) << "Solving Problem:\n" << '\n';
00362 ffout(2) << "\tUnknowns: " << *__unknownList << '\n';
00363 ffout(2) << "\tProblem: " << *__problemExpression << '\n';
00364
00365 if (__domain != 0) {
00366 SolverDriver sd(P, u, discretization,__mesh->mesh(), degreeOfFreedomSet,
00367 SolverDriver::fictitiousFEM);
00368 sd.run();
00369 } else {
00370 SolverDriver sd(P, u, discretization,__mesh->mesh(), degreeOfFreedomSet,
00371 SolverDriver::fem);
00372 sd.run();
00373 }
00374
00375 u.getUserFunction(UL);
00376
00377 UnknownListExpressionSet& UList
00378 = dynamic_cast<UnknownListExpressionSet&>(*__unknownList);
00379 size_t i=0;
00380 for(UnknownListExpressionSet::listType::iterator iunknown = UList.giveList().begin();
00381 iunknown != UList.giveList().end(); ++iunknown) {
00382 UnknownExpression& unknown = **iunknown;
00383
00384 ReferenceCounting<FunctionVariable> unknownVariable
00385 = VariableRepository::instance().findVariable<FunctionVariable>(unknown.name());
00386
00387 const FEMFunctionBase* value = UL[i];
00388 (*unknownVariable)
00389 = new FunctionExpressionValue(value, false);
00390
00391 ++i;
00392 }
00393 }
00394
00395 void SolverExpression::
00396 __solveLegendre(const DiscretizationType& discretization)
00397 {
00398 ReferenceCounting<DegreeOfFreedomSetBuilder> dofBuilder;
00399 if (__domain != 0) {
00400 dofBuilder = new DegreeOfFreedomSetBuilder(discretization,
00401 *__mesh->mesh(),
00402 *__domain->domain());
00403 } else {
00404 dofBuilder = new DegreeOfFreedomSetBuilder(discretization,
00405 *__mesh->mesh());
00406 }
00407
00408 switch (__mesh->mesh()->type()) {
00409 case Mesh::spectralMesh:
00410 case Mesh::octreeMesh: {
00411 break;
00412 }
00413 default: {
00414 throw ErrorHandler(__FILE__,__LINE__,
00415 "cannot use spectral method on '"+__mesh->mesh()->typeName()+"'",
00416 ErrorHandler::normal);
00417 }
00418 }
00419
00420 const DegreeOfFreedomSet& degreeOfFreedomSet
00421 = dofBuilder->degreeOfFreedomSet();
00422
00423 LegendreSolution u(degreeOfFreedomSet);
00424
00425 typedef
00426 std::vector<ReferenceCounting<SpectralFunction> >
00427 UnknownList;
00428 UnknownList UL;
00429
00430 for (UnknownListExpression::iterator i= __unknownList->begin();
00431 i != __unknownList->end(); ++i) {
00432 UnknownExpression& unknown = **i;
00433
00434 ReferenceCounting<FunctionVariable> unknownVariable
00435 = VariableRepository::instance().findVariable<FunctionVariable>(unknown.name());
00436
00437 ConstReferenceCounting<ScalarFunctionBase> f = unknownVariable->expression()->function();
00438
00439 if (f->type() != ScalarFunctionBase::spectral) {
00440 throw ErrorHandler(__FILE__,__LINE__,
00441 "Cannot use non spectral functions as input for spectral method",
00442 ErrorHandler::unexpected);
00443 }
00444 const SpectralFunction* originalSFunction
00445 = dynamic_cast<const SpectralFunction*>(static_cast<const ScalarFunctionBase*>(f));
00446
00447 ReferenceCounting<SpectralFunction> sFunction = new SpectralFunction(*originalSFunction);
00448
00449 UL.push_back(sFunction);
00450 }
00451
00452 u.setUserFunction(UL);
00453
00454 ConstReferenceCounting<Problem> P;
00455 switch (__problemExpression->problemType()) {
00456 case ProblemExpression::pdeSystem: {
00457 P = (PDESystem*)(dynamic_cast<PDESystemExpression&>(*__problemExpression).pdeSystem());
00458 break;
00459 }
00460 case ProblemExpression::variationalProblem: {
00461 P = (VariationalProblem*)(dynamic_cast<VariationalProblemExpression&>(*__problemExpression).variationalProblem());
00462 break;
00463 }
00464 default: {
00465 throw ErrorHandler(__FILE__,__LINE__,
00466 "unexpected problem type",
00467 ErrorHandler::unexpected);
00468 }
00469 }
00470
00471 ffout(2) << "Solving Problem:\n" << '\n';
00472 ffout(2) << "\tUnknowns: " << *__unknownList << '\n';
00473 ffout(2) << "\tProblem: " << *__problemExpression << '\n';
00474
00475 if (__domain != 0) {
00476 SolverDriver sd(P, u, discretization,__mesh->mesh(), degreeOfFreedomSet,
00477 SolverDriver::fictitiousSpectral);
00478 sd.run();
00479 } else {
00480 SolverDriver sd(P, u, discretization,__mesh->mesh(), degreeOfFreedomSet,
00481 SolverDriver::spectral);
00482 sd.run();
00483 }
00484
00485 u.getUserFunction(UL);
00486
00487 UnknownListExpressionSet& UList
00488 = dynamic_cast<UnknownListExpressionSet&>(*__unknownList);
00489 size_t i=0;
00490 for(UnknownListExpressionSet::listType::iterator iunknown = UList.giveList().begin();
00491 iunknown != UList.giveList().end(); ++iunknown) {
00492 UnknownExpression& unknown = **iunknown;
00493
00494 ReferenceCounting<FunctionVariable> unknownVariable
00495 = VariableRepository::instance().findVariable<FunctionVariable>(unknown.name());
00496
00497 const SpectralFunction* value = UL[i];
00498 (*unknownVariable)
00499 = new FunctionExpressionValue(value, false);
00500
00501 ++i;
00502 }
00503 }
