/[ff3d]/ff3d/solver/EliminatedFictitiousDomain.cpp
ViewVC logotype

Diff of /ff3d/solver/EliminatedFictitiousDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by delpinux, Tue Sep 23 19:03:23 2003 UTC revision 1.2 by delpinux, Mon Nov 3 18:56:05 2003 UTC
# Line 17  Line 17 
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());  
 }  
   

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26