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 <[EMAIL PROTECTED]> wrote:
> 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
