Of course putting all my thoughts together in an email made me realize
what I needed to do... all I need to do is use the usual side
assembly... and call elem->is_node_on_side(i,side)... and if it is
then that dof is on the side... then I can just set
Re(i)=soln(dof_indices[i])-bc_value like I want to... overriding
whatever values were assembled into Re before I ever add it into the
residual...
This causes a little bit more work... but not much (especially because
it gets to reuse a bunch of things).
Any other ideas.
Derek
On Thu, Jun 12, 2008 at 1:45 PM, Derek Gaston <[EMAIL PROTECTED]> 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?
>
> 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 u-BC 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(dim-1, 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)=u-bc_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
>
-------------------------------------------------------------------------
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
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
