17 |
|
|
18 |
// $Id$ |
// $Id$ |
19 |
|
|
|
|
|
|
#include <PDESolution.hpp> |
|
|
|
|
20 |
#include <EliminatedFictitiousDomain.hpp> |
#include <EliminatedFictitiousDomain.hpp> |
|
#include <PDEProblem.hpp> |
|
|
|
|
|
#include <Structured3DMesh.hpp> |
|
|
#include <MeshOfHexahedra.hpp> |
|
|
|
|
|
#include <FEMDiscretization.hpp> |
|
|
|
|
21 |
#include <BoundaryConditionDiscretizationElimination.hpp> |
#include <BoundaryConditionDiscretizationElimination.hpp> |
22 |
|
|
|
#include <map> |
|
|
#include <list> |
|
|
|
|
|
#include <KrylovSolver.hpp> |
|
|
|
|
|
#include <MatrixManagement.hpp> |
|
|
|
|
|
#include <SparseMatrix.hpp> |
|
|
|
|
|
#include <Timer.hpp> |
|
|
|
|
23 |
template <typename MeshType> |
template <typename MeshType> |
24 |
void EliminatedFictitiousDomain::__discretize() |
ReferenceCounting<BoundaryConditionDiscretization> |
25 |
|
EliminatedFictitiousDomain::__discretizeBoundaryConditions() |
26 |
{ |
{ |
27 |
MemoryManager MM; |
return new BoundaryConditionDiscretizationElimination<MeshType>(problem(), |
28 |
|
dynamic_cast<MeshType&>(mesh()), |
29 |
bool performAssembling =MM.ReserveMatrix(__A, |
(*__degreeOfFreedomSet)); |
|
problem().numberOfUnknown(), |
|
|
__degreeOfFreedomSet.size()); |
|
|
ffout(2) << "reserve \n"; |
|
|
|
|
|
MM.ReserveVector(__b, |
|
|
problem().numberOfUnknown(), |
|
|
__degreeOfFreedomSet.size()); |
|
|
|
|
|
ffout(2) << "Discretizing the PDE Problem ... "; |
|
|
|
|
|
ReferenceCounting<FEMDiscretization<MeshType> > FEM |
|
|
= new FEMDiscretization<MeshType>(problem(), |
|
|
dynamic_cast<MeshType&>(mesh()), |
|
|
*__A,*__b, __degreeOfFreedomSet); |
|
|
|
|
|
// mesh().CreateReferences(Omega); |
|
|
if (performAssembling) { |
|
|
(*FEM).assembleMatrix(); |
|
|
} else { |
|
|
ffout(2) << "keeping previous operator discretization\n"; |
|
|
} |
|
|
|
|
|
(*FEM).assembleSecondMember(); |
|
|
|
|
|
ffout(2) << "done\n"; |
|
|
|
|
|
ffout(2) << "Discretizing Boundary Conditions\n"; |
|
|
|
|
|
ReferenceCounting<BoundaryConditionDiscretization> bcDiscretization |
|
|
= new BoundaryConditionDiscretizationElimination<MeshType>(problem(), |
|
|
dynamic_cast<MeshType&>(mesh()), |
|
|
__degreeOfFreedomSet); |
|
|
|
|
|
ffout(2) << "Second Member Modification\n"; |
|
|
|
|
|
(*bcDiscretization).setSecondMember(__A,__b); |
|
|
|
|
|
|
|
|
ffout(2) << "Matrix Modification\n"; |
|
|
|
|
|
(*bcDiscretization).setMatrix(__A,__b); |
|
|
|
|
|
ffout(2) << "done\n"; |
|
|
|
|
|
if ((*__A).type() == BaseMatrix::doubleHashedMatrix) { |
|
|
Timer t; |
|
|
t.start(); |
|
|
|
|
|
SparseMatrix* aa |
|
|
= new SparseMatrix(__degreeOfFreedomSet.size(), |
|
|
__degreeOfFreedomSet.numberOfVariables()*27); |
|
|
*aa = static_cast<DoubleHashedMatrix&>(*__A); |
|
|
__A = aa; |
|
|
|
|
|
t.stop(); |
|
|
ffout(2) << "Matrix copy: " << t << '\n'; |
|
|
} |
|
30 |
} |
} |
31 |
|
|
32 |
void EliminatedFictitiousDomain::Discretize (ConstReferenceCounting<Problem> givenProblem) |
ReferenceCounting<BoundaryConditionDiscretization> |
33 |
|
EliminatedFictitiousDomain::discretizeBoundaryConditions() |
34 |
{ |
{ |
|
ffout(2) << "Adding Characteristic function of the domain to PDEs ... " << std::flush; |
|
|
ConstReferenceCounting<UserFunction> u = new OneDomainUserFunction((*givenProblem).domain()); |
|
|
__problem = (*givenProblem) * u; |
|
|
ffout(2) << "done\n"; |
|
|
|
|
35 |
switch (mesh().type()) { |
switch (mesh().type()) { |
36 |
case Mesh::cartesianHexahedraMesh: { |
case Mesh::cartesianHexahedraMesh: { |
37 |
this->__discretize<Structured3DMesh>(); |
return this->__discretizeBoundaryConditions<Structured3DMesh>(); |
38 |
break; |
break; |
39 |
} |
} |
40 |
default: { |
default: { |
41 |
fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n"; |
fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n"; |
42 |
fferr(0) << "this mesh type is not supported by FDM elimination\n"; |
fferr(0) << "this mesh type is not supported by FDM penalty\n"; |
43 |
std::exit(1); |
std::exit(1); |
44 |
|
return 0; |
45 |
} |
} |
46 |
} |
} |
47 |
} |
} |
|
|
|
|
void EliminatedFictitiousDomain::Compute (Solution& U) |
|
|
{ |
|
|
PDESolution& u = static_cast<PDESolution&>(U); |
|
|
|
|
|
ffout(2) << "Solving the linear system\n"; |
|
|
#warning BAD BAD BAD ... |
|
|
|
|
|
KrylovSolver K(*__A, *__b, __degreeOfFreedomSet, |
|
|
new Structured3DMeshShape(dynamic_cast<Structured3DMesh&>(mesh()).shape())); |
|
|
|
|
|
K.solve(u.values(), problem()); |
|
|
} |
|
|
|
|