// This file is part of ff3d - http://www.freefem.org/ff3d // Copyright (C) 2001, 2002, 2003 Stéphane Del Pino // This program is free software; you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation; either version 2, or (at your option) // any later version. // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software Foundation, // Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // $Id: EliminatedFictitiousDomain.cpp,v 1.1 2003/09/23 19:03:23 delpinux Exp $ #include #include #include #include #include #include #include #include #include #include #include #include #include template void EliminatedFictitiousDomain::__discretize() { MemoryManager MM; bool performAssembling =MM.ReserveMatrix(__A, problem().numberOfUnknown(), __degreeOfFreedomSet.size()); ffout(2) << "reserve \n"; MM.ReserveVector(__b, problem().numberOfUnknown(), __degreeOfFreedomSet.size()); ffout(2) << "Discretizing the PDE Problem ... "; ReferenceCounting > FEM = new FEMDiscretization(problem(), dynamic_cast(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 bcDiscretization = new BoundaryConditionDiscretizationElimination(problem(), dynamic_cast(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(*__A); __A = aa; t.stop(); ffout(2) << "Matrix copy: " << t << '\n'; } } void EliminatedFictitiousDomain::Discretize (ConstReferenceCounting givenProblem) { ffout(2) << "Adding Characteristic function of the domain to PDEs ... " << std::flush; ConstReferenceCounting u = new OneDomainUserFunction((*givenProblem).domain()); __problem = (*givenProblem) * u; ffout(2) << "done\n"; switch (mesh().type()) { case Mesh::cartesianHexahedraMesh: { this->__discretize(); break; } default: { fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n"; fferr(0) << "this mesh type is not supported by FDM elimination\n"; std::exit(1); } } } void EliminatedFictitiousDomain::Compute (Solution& U) { PDESolution& u = static_cast(U); ffout(2) << "Solving the linear system\n"; #warning BAD BAD BAD ... KrylovSolver K(*__A, *__b, __degreeOfFreedomSet, new Structured3DMeshShape(dynamic_cast(mesh()).shape())); K.solve(u.values(), problem()); }