From: Derek Gaston <friedmud@gm...>  20080612 19:45:52

I'm in an interesting position where I need to do a numeric_vector.set() after doing a bunch of numeric_vector.add() operations... and of course PETSc doesn't like that (says that the object is in the wrong state). Doing the add and set in serial works fine... What are my options here? Alternatively, I might be going about the whole thing incorrectly. The deal is that I need to modify my residual in my jacobian free method to account for Dirichlet B.Cs and I want to do it strongly (ie put uBC for the residual for those dofs). I can't seem to figure out how to do in the "normal" libMesh assembly way of doing things. We usually integrate these B.C.s over the entire element on the edge and take advantage of the fact that the support for the interior dofs is usually zero on that edge. This means that I can't use a loop like this: ###### AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qface(dim1, FIFTH); ... for ( ; el != end_el ; ++el) { ... for (unsigned int qp=0; qp<qrule.n_points(); qp++) { //Interior Residual assembly } for (unsigned int side=0; side<elem>n_sides(); side++) if (elem>neighbor(side) == NULL) for (unsigned int qp=0; qp<qface.n_points(); qp++) for (unsigned int i=0; i<phi_face.size(); i++) Re(i)=ubc_value; ###### Because the Re(i)'s include interior degrees of freedom. What I've done to get around this feels all kinds of hackish... it goes something like this: ###### // Do all interior residual calculations... residual.add_vector (Re, dof_indices); for (unsigned int side=0; side<elem>n_sides(); side++) if (elem>neighbor(side) == NULL) { unsigned int boundary_id = mesh.boundary_info>boundary_id (elem, side); std::vector<unsigned int> side_dof_indices; AutoPtr<Elem> side_elem = elem>build_side(side); dof_map.dof_indices (side_elem.get(), side_dof_indices); side_Re.resize(side_dof_indices.size()); for(unsigned int j=0; j<side_dof_indices.size(); j++) side_Re(j)=soln(side_dof_indices[j])bc_value; } residual.insert(side_Re, side_dof_indices); ###### (Note, the above is missing a bunch of pieces that are actually in my code... this is just to give you a flavor) This works great in serial... but not in parallel where PETSc doesn't like the combo of add_vector() and insert().... Is there a better way to do this? Is there a way to tell if the dofs are on the boundary or not without having to get the side_dof_indices? That would allow me to modify Re before I ever add it to the residual... Thanks for any help! Derek 