From: Derek G. <fri...@gm...> - 2008-06-12 19:45:52
|
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 |
From: John P. <jwp...@gm...> - 2008-06-12 20:00:04
|
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 > |
From: Derek G. <fri...@gm...> - 2008-06-12 20:02:35
|
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 >> > |
From: Derek G. <fri...@gm...> - 2008-06-12 20:01:13
|
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 <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? > > 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 > |
From: Roy S. <ro...@st...> - 2008-06-12 20:08:34
|
On Thu, 12 Jun 2008, Derek Gaston wrote: > 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 looks good, but you might find it easier in the long run to use FE::dofs_on_side() rather than Elem::is_node_on_side; the latter might make it too easy to write code that only works on isoparametric Lagrange elements. --- Roy |
From: Derek G. <fri...@gm...> - 2008-06-12 20:12:58
|
It's kind of hard to see how to go from dofs_on_side back to the mapping to Re()... you have to know what "node number" those dofs map to for the current element. Or am I missing something? Derek On Thu, Jun 12, 2008 at 2:08 PM, Roy Stogner <ro...@st...> wrote: > > On Thu, 12 Jun 2008, Derek Gaston wrote: > >> 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 looks good, but you might find it easier in the long run to use > FE::dofs_on_side() rather than Elem::is_node_on_side; the latter might > make it too easy to write code that only works on isoparametric > Lagrange elements. > --- > Roy > |
From: Benjamin K. <ben...@na...> - 2008-06-12 20:25:04
|
Many, many ages ago we did Lagrange BCs ('cus there weren't any other types of elements yet!) strongly. The trick is to do it at the element matrix/vector level before inserting into the global matrix. There is actually a member function in DenseMatrix to do this: Ke.condense(i,j,val,Fe); is what you want. This will impose val strongly by doing the right thing at the element level. Of course you can add these and instead of getting val=foo at the linear system level you get N*val = N*foo There is another, important issue though that John learned via suboptimal convergence on triangular meshes... Usually we only concern ourselves with BCs on elements whose sides intersect the boundary, e.g. NULL neighbors... If you are going to impose the value strongly through the element matrix approach you have to do it for all elements which touch the DOF. Finally, you can effectively accomplish the same thing with penalty using add() instead of set(), just add 1e30 to the global matrix diagonal and 1e30*value to the RHS, all the other entries will "go away" in floating point. -Ben On 6/12/08 3:12 PM, "Derek Gaston" <fri...@gm...> wrote: > It's kind of hard to see how to go from dofs_on_side back to the > mapping to Re()... you have to know what "node number" those dofs map > to for the current element. Or am I missing something? > > Derek > > On Thu, Jun 12, 2008 at 2:08 PM, Roy Stogner <ro...@st...> wrote: >> >> On Thu, 12 Jun 2008, Derek Gaston wrote: >> >>> 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 looks good, but you might find it easier in the long run to use >> FE::dofs_on_side() rather than Elem::is_node_on_side; the latter might >> make it too easy to write code that only works on isoparametric >> Lagrange elements. >> --- >> Roy >> > > ------------------------------------------------------------------------- > 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 |
From: Derek G. <fri...@gm...> - 2008-06-12 20:46:39
|
Note that we don't have a matrix to condense in our matrix free scheme ;-) I ultimately think that we're going to use a BoundaryMesh to do this... we'll just loop over the nodes in the boundary mesh and set their residual as a post-processing step.... and avoid all of these problems (including the triangle/tet problem). Derek On Thu, Jun 12, 2008 at 2:24 PM, Benjamin Kirk <ben...@na...> wrote: > Many, many ages ago we did Lagrange BCs ('cus there weren't any other types > of elements yet!) strongly. The trick is to do it at the element > matrix/vector level before inserting into the global matrix. > > There is actually a member function in DenseMatrix to do this: > > Ke.condense(i,j,val,Fe); > > is what you want. This will impose val strongly by doing the right thing at > the element level. Of course you can add these and instead of getting > > val=foo > > at the linear system level you get > > N*val = N*foo > > > There is another, important issue though that John learned via suboptimal > convergence on triangular meshes... Usually we only concern ourselves with > BCs on elements whose sides intersect the boundary, e.g. NULL neighbors... > If you are going to impose the value strongly through the element matrix > approach you have to do it for all elements which touch the DOF. > > > Finally, you can effectively accomplish the same thing with penalty using > add() instead of set(), just add 1e30 to the global matrix diagonal and > 1e30*value to the RHS, all the other entries will "go away" in floating > point. > > -Ben > > > > On 6/12/08 3:12 PM, "Derek Gaston" <fri...@gm...> wrote: > >> It's kind of hard to see how to go from dofs_on_side back to the >> mapping to Re()... you have to know what "node number" those dofs map >> to for the current element. Or am I missing something? >> >> Derek >> >> On Thu, Jun 12, 2008 at 2:08 PM, Roy Stogner <ro...@st...> wrote: >>> >>> On Thu, 12 Jun 2008, Derek Gaston wrote: >>> >>>> 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 looks good, but you might find it easier in the long run to use >>> FE::dofs_on_side() rather than Elem::is_node_on_side; the latter might >>> make it too easy to write code that only works on isoparametric >>> Lagrange elements. >>> --- >>> Roy >>> >> >> ------------------------------------------------------------------------- >> 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 > > |
From: John P. <jwp...@gm...> - 2008-06-12 21:02:14
|
On Thu, Jun 12, 2008 at 3:46 PM, Derek Gaston <fri...@gm...> wrote: > Note that we don't have a matrix to condense in our matrix free scheme ;-) > > I ultimately think that we're going to use a BoundaryMesh to do > this... we'll just loop over the nodes in the boundary mesh and set > their residual as a post-processing step.... and avoid all of these > problems (including the triangle/tet problem). I know this is still in the idea stage, but I thought I would mention -- currently when you use the BoundaryInfo object to sync() with a BoundaryMesh, the BoundaryMesh makes its own *copy* of the nodes from the original mesh. (The reason for this was so that the BoundaryMesh could survive on its own, even if the original Mesh was destroyed. Also you can alter the BoundaryMesh independently of the original if you desire.) A better approach might be to define a boundary node_iterator in the style of the other node iterators... but I'm sure you will cross this bridge when you come to it. -J > > Derek > > On Thu, Jun 12, 2008 at 2:24 PM, Benjamin Kirk <ben...@na...> wrote: >> Many, many ages ago we did Lagrange BCs ('cus there weren't any other types >> of elements yet!) strongly. The trick is to do it at the element >> matrix/vector level before inserting into the global matrix. >> >> There is actually a member function in DenseMatrix to do this: >> >> Ke.condense(i,j,val,Fe); >> >> is what you want. This will impose val strongly by doing the right thing at >> the element level. Of course you can add these and instead of getting >> >> val=foo >> >> at the linear system level you get >> >> N*val = N*foo >> >> >> There is another, important issue though that John learned via suboptimal >> convergence on triangular meshes... Usually we only concern ourselves with >> BCs on elements whose sides intersect the boundary, e.g. NULL neighbors... >> If you are going to impose the value strongly through the element matrix >> approach you have to do it for all elements which touch the DOF. >> >> >> Finally, you can effectively accomplish the same thing with penalty using >> add() instead of set(), just add 1e30 to the global matrix diagonal and >> 1e30*value to the RHS, all the other entries will "go away" in floating >> point. >> >> -Ben >> >> >> >> On 6/12/08 3:12 PM, "Derek Gaston" <fri...@gm...> wrote: >> >>> It's kind of hard to see how to go from dofs_on_side back to the >>> mapping to Re()... you have to know what "node number" those dofs map >>> to for the current element. Or am I missing something? >>> >>> Derek >>> >>> On Thu, Jun 12, 2008 at 2:08 PM, Roy Stogner <ro...@st...> wrote: >>>> >>>> On Thu, 12 Jun 2008, Derek Gaston wrote: >>>> >>>>> 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 looks good, but you might find it easier in the long run to use >>>> FE::dofs_on_side() rather than Elem::is_node_on_side; the latter might >>>> make it too easy to write code that only works on isoparametric >>>> Lagrange elements. >>>> --- >>>> Roy >>>> >>> >>> ------------------------------------------------------------------------- >>> 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 >> >> > > ------------------------------------------------------------------------- > 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 > |
From: Derek G. <fri...@gm...> - 2008-06-12 21:19:11
|
Ah - that is good info... I thought the BoundaryMesh just had pointers to the other mesh's data structures... kind of blows that idea out of the water... So how are other people doing boundary conditions with tri's and tets? With Dirichlet you can just use a penalty to swamp everything out. But with Nuemann? Derek On Jun 12, 2008, at 3:02 PM, John Peterson wrote: > On Thu, Jun 12, 2008 at 3:46 PM, Derek Gaston <fri...@gm...> > wrote: >> Note that we don't have a matrix to condense in our matrix free >> scheme ;-) >> >> I ultimately think that we're going to use a BoundaryMesh to do >> this... we'll just loop over the nodes in the boundary mesh and set >> their residual as a post-processing step.... and avoid all of these >> problems (including the triangle/tet problem). > > I know this is still in the idea stage, but I thought I would mention > -- currently when you use the BoundaryInfo object to sync() with a > BoundaryMesh, the BoundaryMesh makes its own *copy* of the nodes from > the original mesh. (The reason for this was so that the BoundaryMesh > could survive on its own, even if the original Mesh was destroyed. > Also you can alter the BoundaryMesh independently of the original if > you desire.) > > A better approach might be to define a boundary node_iterator in the > style of the other node iterators... but I'm sure you will cross this > bridge when you come to it. > > -J > > > >> >> Derek >> >> On Thu, Jun 12, 2008 at 2:24 PM, Benjamin Kirk <ben...@na... >> > wrote: >>> Many, many ages ago we did Lagrange BCs ('cus there weren't any >>> other types >>> of elements yet!) strongly. The trick is to do it at the element >>> matrix/vector level before inserting into the global matrix. >>> >>> There is actually a member function in DenseMatrix to do this: >>> >>> Ke.condense(i,j,val,Fe); >>> >>> is what you want. This will impose val strongly by doing the >>> right thing at >>> the element level. Of course you can add these and instead of >>> getting >>> >>> val=foo >>> >>> at the linear system level you get >>> >>> N*val = N*foo >>> >>> >>> There is another, important issue though that John learned via >>> suboptimal >>> convergence on triangular meshes... Usually we only concern >>> ourselves with >>> BCs on elements whose sides intersect the boundary, e.g. NULL >>> neighbors... >>> If you are going to impose the value strongly through the element >>> matrix >>> approach you have to do it for all elements which touch the DOF. >>> >>> >>> Finally, you can effectively accomplish the same thing with >>> penalty using >>> add() instead of set(), just add 1e30 to the global matrix >>> diagonal and >>> 1e30*value to the RHS, all the other entries will "go away" in >>> floating >>> point. >>> >>> -Ben >>> >>> >>> >>> On 6/12/08 3:12 PM, "Derek Gaston" <fri...@gm...> wrote: >>> >>>> It's kind of hard to see how to go from dofs_on_side back to the >>>> mapping to Re()... you have to know what "node number" those dofs >>>> map >>>> to for the current element. Or am I missing something? >>>> >>>> Derek >>>> >>>> On Thu, Jun 12, 2008 at 2:08 PM, Roy Stogner <ro...@st...> >>>> wrote: >>>>> >>>>> On Thu, 12 Jun 2008, Derek Gaston wrote: >>>>> >>>>>> 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 looks good, but you might find it easier in the long run to >>>>> use >>>>> FE::dofs_on_side() rather than Elem::is_node_on_side; the latter >>>>> might >>>>> make it too easy to write code that only works on isoparametric >>>>> Lagrange elements. >>>>> --- >>>>> Roy >>>>> >>>> >>>> ------------------------------------------------------------------------- >>>> 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 >>> >>> >> >> ------------------------------------------------------------------------- >> 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 >> |
From: John P. <jwp...@gm...> - 2008-06-12 21:25:22
|
On Thu, Jun 12, 2008 at 4:18 PM, Derek Gaston <fri...@gm...> wrote: > Ah - that is good info... I thought the BoundaryMesh just had pointers to > the other mesh's data structures... kind of blows that idea out of the > water... > > So how are other people doing boundary conditions with tri's and tets? With > Dirichlet you can just use a penalty to swamp everything out. But with > Nuemann? Neumann bcs are already in terms of an edge (or face) integral. So these elements with zero support (only a node) on the boundary don't figure into the Neumann case anyway. -J > Derek > > On Jun 12, 2008, at 3:02 PM, John Peterson wrote: > >> On Thu, Jun 12, 2008 at 3:46 PM, Derek Gaston <fri...@gm...> wrote: >>> >>> Note that we don't have a matrix to condense in our matrix free scheme >>> ;-) >>> >>> I ultimately think that we're going to use a BoundaryMesh to do >>> this... we'll just loop over the nodes in the boundary mesh and set >>> their residual as a post-processing step.... and avoid all of these >>> problems (including the triangle/tet problem). >> >> I know this is still in the idea stage, but I thought I would mention >> -- currently when you use the BoundaryInfo object to sync() with a >> BoundaryMesh, the BoundaryMesh makes its own *copy* of the nodes from >> the original mesh. (The reason for this was so that the BoundaryMesh >> could survive on its own, even if the original Mesh was destroyed. >> Also you can alter the BoundaryMesh independently of the original if >> you desire.) >> >> A better approach might be to define a boundary node_iterator in the >> style of the other node iterators... but I'm sure you will cross this >> bridge when you come to it. >> >> -J >> >> >> >>> >>> Derek >>> >>> On Thu, Jun 12, 2008 at 2:24 PM, Benjamin Kirk <ben...@na...> >>> wrote: >>>> >>>> Many, many ages ago we did Lagrange BCs ('cus there weren't any other >>>> types >>>> of elements yet!) strongly. The trick is to do it at the element >>>> matrix/vector level before inserting into the global matrix. >>>> >>>> There is actually a member function in DenseMatrix to do this: >>>> >>>> Ke.condense(i,j,val,Fe); >>>> >>>> is what you want. This will impose val strongly by doing the right >>>> thing at >>>> the element level. Of course you can add these and instead of getting >>>> >>>> val=foo >>>> >>>> at the linear system level you get >>>> >>>> N*val = N*foo >>>> >>>> >>>> There is another, important issue though that John learned via >>>> suboptimal >>>> convergence on triangular meshes... Usually we only concern ourselves >>>> with >>>> BCs on elements whose sides intersect the boundary, e.g. NULL >>>> neighbors... >>>> If you are going to impose the value strongly through the element matrix >>>> approach you have to do it for all elements which touch the DOF. >>>> >>>> >>>> Finally, you can effectively accomplish the same thing with penalty >>>> using >>>> add() instead of set(), just add 1e30 to the global matrix diagonal and >>>> 1e30*value to the RHS, all the other entries will "go away" in floating >>>> point. >>>> >>>> -Ben >>>> >>>> >>>> >>>> On 6/12/08 3:12 PM, "Derek Gaston" <fri...@gm...> wrote: >>>> >>>>> It's kind of hard to see how to go from dofs_on_side back to the >>>>> mapping to Re()... you have to know what "node number" those dofs map >>>>> to for the current element. Or am I missing something? >>>>> >>>>> Derek >>>>> >>>>> On Thu, Jun 12, 2008 at 2:08 PM, Roy Stogner <ro...@st...> wrote: >>>>>> >>>>>> On Thu, 12 Jun 2008, Derek Gaston wrote: >>>>>> >>>>>>> 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 looks good, but you might find it easier in the long run to use >>>>>> FE::dofs_on_side() rather than Elem::is_node_on_side; the latter might >>>>>> make it too easy to write code that only works on isoparametric >>>>>> Lagrange elements. >>>>>> --- >>>>>> Roy >>>>>> >>>>> >>>>> >>>>> ------------------------------------------------------------------------- >>>>> 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 >>>> >>>> >>> >>> ------------------------------------------------------------------------- >>> 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 >>> > > |
From: Benjamin K. <ben...@na...> - 2008-06-12 21:25:57
|
> So how are other people doing boundary conditions with tri's and > tets? With Dirichlet you can just use a penalty to swamp everything > out. But with Nuemann? Neumann is not an issue because the BC is in terms of the boundary integral, so you only consider it with elements whose sides intersect the boundary. The elements who only touch the boundary through a node do not enter into the Neumann BC. |
From: Derek G. <fri...@gm...> - 2008-06-12 21:33:34
|
Right - maybe I should think about these things before I say them ;-) So now the question is... if you were in my place and wanted to set a u-bc_val type residual for boundary elements... what would you do for tri's and tets? It seems like a post-processing step is the best way (so that something from an element with only a node on the boundary doesn't have a chance to contribute anything). What would that post-processing step be? Are there any facilities already in libMesh that could help out? I mean, I can certainly loop over the local_active_elements and add all the nodes on the boundaries to a std::vector<Node *> right after I read the mesh... but there is a bit of a problem with adaptivity in this case... and it feels hackish. But I suppose I was hoping that this is what a BoundaryMesh was... so I guess it's not all that much different. Thanks, Derek On Thu, Jun 12, 2008 at 3:25 PM, Benjamin Kirk <ben...@na...> wrote: >> So how are other people doing boundary conditions with tri's and >> tets? With Dirichlet you can just use a penalty to swamp everything >> out. But with Nuemann? > > Neumann is not an issue because the BC is in terms of the boundary integral, > so you only consider it with elements whose sides intersect the boundary. > The elements who only touch the boundary through a node do not enter into > the Neumann BC. > > |
From: John P. <jwp...@gm...> - 2008-06-12 21:48:15
|
On Thu, Jun 12, 2008 at 4:33 PM, Derek Gaston <fri...@gm...> wrote: > Right - maybe I should think about these things before I say them ;-) > > So now the question is... if you were in my place and wanted to set a > u-bc_val type residual for boundary elements... what would you do for > tri's and tets? It seems like a post-processing step is the best way > (so that something from an element with only a node on the boundary > doesn't have a chance to contribute anything). What would that > post-processing step be? Are there any facilities already in libMesh > that could help out? > > I mean, I can certainly loop over the local_active_elements and add > all the nodes on the boundaries to a std::vector<Node *> right after I > read the mesh... but there is a bit of a problem with adaptivity in > this case... and it feels hackish. But I suppose I was hoping that > this is what a BoundaryMesh was... so I guess it's not all that much > different. BoundaryInfo::build_node_list() does what you suggest. Given the current capabilities of the library, I'd say rebuilding this list after each assembly, and going through and modifying those entries in the residual vector is the way to do it. This should also be done in parallel, i.e. you only set the dof indexes from these nodes which are on your processor. -J > > Thanks, > Derek > > On Thu, Jun 12, 2008 at 3:25 PM, Benjamin Kirk <ben...@na...> wrote: >>> So how are other people doing boundary conditions with tri's and >>> tets? With Dirichlet you can just use a penalty to swamp everything >>> out. But with Nuemann? >> >> Neumann is not an issue because the BC is in terms of the boundary integral, >> so you only consider it with elements whose sides intersect the boundary. >> The elements who only touch the boundary through a node do not enter into >> the Neumann BC. >> >> > |
From: Derek G. <fri...@gm...> - 2008-06-12 21:49:06
|
Great - that's exactly what I needed! Thanks John! Derek On Thu, Jun 12, 2008 at 3:48 PM, John Peterson <jwp...@gm...> wrote: > On Thu, Jun 12, 2008 at 4:33 PM, Derek Gaston <fri...@gm...> wrote: >> Right - maybe I should think about these things before I say them ;-) >> >> So now the question is... if you were in my place and wanted to set a >> u-bc_val type residual for boundary elements... what would you do for >> tri's and tets? It seems like a post-processing step is the best way >> (so that something from an element with only a node on the boundary >> doesn't have a chance to contribute anything). What would that >> post-processing step be? Are there any facilities already in libMesh >> that could help out? >> >> I mean, I can certainly loop over the local_active_elements and add >> all the nodes on the boundaries to a std::vector<Node *> right after I >> read the mesh... but there is a bit of a problem with adaptivity in >> this case... and it feels hackish. But I suppose I was hoping that >> this is what a BoundaryMesh was... so I guess it's not all that much >> different. > > BoundaryInfo::build_node_list() does what you suggest. Given the > current capabilities of the library, I'd say rebuilding this list > after each assembly, and going through and modifying those entries in > the residual vector is the way to do it. This should also be done in > parallel, i.e. you only set the dof indexes from these nodes which are > on your processor. > > -J > > > > >> >> Thanks, >> Derek >> >> On Thu, Jun 12, 2008 at 3:25 PM, Benjamin Kirk <ben...@na...> wrote: >>>> So how are other people doing boundary conditions with tri's and >>>> tets? With Dirichlet you can just use a penalty to swamp everything >>>> out. But with Nuemann? >>> >>> Neumann is not an issue because the BC is in terms of the boundary integral, >>> so you only consider it with elements whose sides intersect the boundary. >>> The elements who only touch the boundary through a node do not enter into >>> the Neumann BC. >>> >>> >> > |
From: Roy S. <ro...@st...> - 2008-06-12 20:33:06
|
On Thu, 12 Jun 2008, Roy Stogner wrote: > On Thu, 12 Jun 2008, Derek Gaston wrote: > >> It's kind of hard to see how to go from dofs_on_side back to the >> mapping to Re()... you have to know what "node number" those dofs map >> to for the current element. Or am I missing something? > > IIRC dofs_on_side gives you local dof numbers; you still need the > result of a dof_indices call to look up the corresponding global dof > indices. Oh, and to really be general this implies that your residual terms aren't just the pointwise difference between current value and desired value at node i, they're the sum of the side integral of the difference between current and desired values weighted by test function i. So basically you still evaluate those integrals, but you use the dofs_on_side information to decide which test functions to skip. --- Roy |