On Thu, Jun 12, 2008 at 2: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?
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 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...
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
> _______________________________________________
> Libmesh-users mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
-------------------------------------------------------------------------
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
