Re: [Libmesh-users] _Set_ing values in numeric_vectors (or doing strong Dirichlet BCs in matrix free mode)

[Derek G. <fri...@gm...>]

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 <jwp...@gm...> wrote:
> On Thu, Jun 12, 2008 at 2:45 PM, Derek Gaston <fri...@gm...> 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
>> Lib...@li...
>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>
>