On Thu, 9 Sep 2010, Tim Kroeger wrote: > On Thu, 9 Sep 2010, Jed Brown wrote: > >> On Thu, 9 Sep 2010 15:27:42 +0200 (CEST), Tim Kroeger >> <[email protected]> wrote: >>> Dear all, >>> >>> Let X be a computational domain, covered by a Mesh on which an >>> EquationSystems object is based on. Let X2 be a subset of X, given by >>> the union of selected grid elements (which are, say, marked by certain >>> values of subdomain_id). >>> >>> Assume I want (at some point in my application) to solve a system only >>> on X2 (say, using Dirichlet boundary conditions on all boundaries of >>> X2 that are not boundaries of X). >>> >>> I can easily achieve this by assembling Dirichlet conditions >>> everywhere ouside X2 and then solving as usual. However, then I >>> cannot benefit from the performance gain that I should theoretically >>> have if X2 contains much less elements than X. This is in particular >>> true if I am using a direct solver (such as SuperLU_DIST, wrapped via >>> PETSc). >> >> Do you want to assemble X, or are you really only working with X2? > > Most probably, my assmble method will loop over X but skip every > element that is not contained in X2. This seems to be the easiest > assemble method at the moment. > >> If the former, you can MatGetSubMatrix (you just need an index set >> defining the subdomain) the part you want and solve with that. > > That sounds like a good idea. I will think about this. It should > certainly be worth implementing a general method for this in libMesh > somehow (with an exception trown if a solver other than PETSc is > used...). The question is what the API should look like.
My suggestion for an API in libMesh is as follows: LinearSolver gets a method LinearSolver::solve_only_on(const std::set<unsigned int>*const) after a call to which each call to LinearSolver::solve() will remove all dofs not contained in the list from the matrix before actually solving. The default implementation of this will be libmesh_not_implemented(), and only for PetscLinearSolver, I will implement the solution using MatGetSubMatrix() as Jed pointed out. (Optionally, an efficient load balancing using PCFieldSplit() can be added later.) Likewise, System gets a method System::solve_only_on(const std::set<subdomain_id_type>*const) which makes System::solve() call LinearSolver::solve_only_on() with the correct index set, that is, sucht that it solves on all dofs that are associated with at least one elem having a subdomain_id value that is contained in the list. Calling either solve_only_on(NULL) should in both cases restore the default behaviour. libMesh developers, what do you think? Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236 ------------------------------------------------------------------------------ This SF.net Dev2Dev email is sponsored by: Show off your parallel programming skills. Enter the Intel(R) Threading Challenge 2010. http://p.sf.net/sfu/intel-thread-sfd _______________________________________________ Libmesh-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
