From: Derek Gaston <friedmud@gm...>  20080612 20:02:35

Thanks for the reply John. I was hesitant to call close()... I wasn't entirely sure that will work correctly with NonlinearImplicitSystem. But I think I've found another way to do it. Derek On Thu, Jun 12, 2008 at 2:00 PM, John Peterson <jwpeterson@...> wrote: > On Thu, Jun 12, 2008 at 2:45 PM, Derek Gaston <friedmud@...> wrote: >> 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? > > Did you try calling close() before trying to set() values? close > basically calls VectorAssemblyBegin() and End() so that afterward, the > object will be in the right "state". > >> 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... > > Well, once you have a pointer to a node that you know is on the > boundary (should be able to get this information in a variety of ways) > you can do: > > mesh.node_ptr(*it)>dof_number(0, /*system #*/ > v, /*var # */ > 0); /*component #*/ > > to find the global degree of freedom you'd need to modify by hand to > impose the Dirichlet BC. This is specific to Lagrange elements. > > > J > > >> >> Thanks for any help! >> >> Derek >> >>  >> Check out the new SourceForge.net Marketplace. >> It's the best place to buy or sell services for >> just about anything Open Source. >> http://sourceforge.net/services/buy/index.php >> _______________________________________________ >> Libmeshusers mailing list >> Libmeshusers@... >> https://lists.sourceforge.net/lists/listinfo/libmeshusers >> > 