From: Tim Kroeger <tim.kroeger@ce...>  20081023 13:12:17

Dear John/Jed/Roy, On Wed, 22 Oct 2008, John Peterson wrote: >> You might later implement shell matrices for other solver >> packages, and I would find it more appropriate to have a DyadMatrix class as >> the user API and let everything that depends on the solver package be hidden >> from the user. > > This is a good point and raises the issue of how one will interact > with these new matrix types in a general System setting. Up to this > point your modifications have been to the PetscLinearSolver but > eventually I would like this to work within the confines of the > existing ImplicitSystem. > > My current thought is to extend ImplicitSystem::add_matrix(string) to > also take a NumericMatrix* (default NULL). Current users can then > continue with the default SparseMatrix::build behavior currently in > add_matrix(), while giving the opportunity for other users to attach > any other kind of ShellMatrix at the same time. One just needs to > remember that the System will finally delete any Matrix added to it. That's a good idea. You should then remember that there should also be some mechanism to make ImplicitSystem::matrix something different than a SparseMatrix. This way, my method LinearImplicitSystem::attach_shell_matrix() is not required any more. > To address the issue of the shell matrices being tied to Petsc, see > the attached new class diagram. I think we should continue using > the pseudovirtualconstructor idea which is now employed in > SparseMatrix::build() which basically gives you the correct type > depending on what solver package you have asked for. All the Petsc > types are now multiply derived from (and dynamic_cast'able to) > PetscMatrixBase, which will be little more than a C++ wrapper for a > Petsc Mat object. There is no need for the multipleinheritance > diamond in this hierarchy. I like your new design much more than the previous one. However, a user that needs a new version of a shell matrix will have to do some PETSc programming now, which wouldn't be required by my implementation. Well, that might be acceptable. What do you think about Jed's idea? >> By the way: What about my question concerning matrix/vector contraining? I >> didn't receive any answer to that from any of you guys. > > I'm not quite sure what the answer is, perhaps someone who understands > the constraints better will eventually comment. I think it must > depend on what the dyad matrix represents, and how it should be > constrained? Since nobody replied, I think I have to find out myself. I see that all the constrain methods build the socalled constraint matrix C (which is not necessarily quadratic), and then: * for a vector v (right hand side), they compute C'*v (where C' denotes the transposed matrix) and clear the components that correspond to contrained dofs; * for a matrix K, they compute C'*K*C and clear the rows (but not the columns) that correspond to contrained dofs except for the diagonal entries, which are set to 1, and except for certain other entries of these rows, which are set to certain other values. Since a DyadMatrix will in no case be used as the system matrix alone (because it's highly singular), it seems appropriate to constrain B=v*w' by computing C'*v and C'*w and then clearing the components of v (but not of w) that correspond to constrained dofs. I now implemented DofMap::constrain_element_dyad_matrix() which does just that (see attached patch), but I'm not yet satisfied with the user interface since it is now necessary to take of copy of the dof_indices inside the assemble function. I wonder what a better user interface would be. Of course, I could add DofMap::constrain_element_matrix_and_vector_and_dyad_matrix(), but what if the next user needs two dyadic matrices or whatever. At the moment, I would like to unify all the constraining methods to something like DofMap::constrain(std::vector<DenseMatrix<Number> >& matrices, std::vector<DenseVector<Number> >& vectors, std::vector<std::pair<DenseVector<Number>,DenseVector<Number> > >& dyad_matrices, std::vector<unsigned int>& dofs, bool symmetric_constraint_rows) and replace the existing methods with the appropriate calls to that function. Well, in order not to copy marices and vectors, it would be more sensible to use vectors of pointers, i.e. DofMap::constrain(std::vector<DenseMatrix<Number>*>& matrices, std::vector<DenseVector<Number>*>& vectors, std::vector<std::pair<DenseVector<Number>*,DenseVector<Number>*> >& dyad_matrices, std::vector<unsigned int>& dofs, bool symmetric_constraint_rows) What do you think? Of course, there should be another function similar to this but allowing for different dofs in rows and columns (like it exists for constrain_element_matrix()). By the way: I will be out of office on friday, so don't be surprised if I don't reply again before monday. Best Regards, Tim  Dr. Tim Kroeger Phone +494212187710 tim.kroeger@..., tim.kroeger@... Fax +494212184236 MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany Amtsgericht Bremen HRB 16222 Geschaeftsfuehrer: Prof. Dr. H.O. Peitgen 