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
