## Re: [Libmesh-users] Question about Dirichlet BC imposed by penaltymethod

 Re: [Libmesh-users] Question about Dirichlet BC imposed by penaltymethod From: Kirk, Benjamin \(JSC-EG\) - 2006-06-23 17:18:40 ```What you want is the linear system to reduce to the L2 projection of the boundary data for the boundary degrees of freedom. The way we tend to do that is to add the L2 projection, scaled by some *huge* penalty value, to the matrix and RHS. In floating point arithmetic the idea is that you add these huge values to the matrix and the other values in the matrix are no longer significant and essentially are 0. One reason for doing this is that you only need to assess boundary conditions in elements that have a full DIM-1 face on the boundary and can ignore elements with lower-dimensional boundary traces (edges in 3D, nodes in 2D/3D). =20 I am not familiar with your formulation, but I am almost certain you want the penalty value multiplying the Me() term and not the Ke() term. Again, the notion is to add the L2 projection of the boundary data, and the mass matrix contains the terms necessary for the L2 projection on the matrix side. -Ben =20 -----Original Message----- From: libmesh-users-bounces@... [mailto:libmesh-users-bounces@...] On Behalf Of David Xu Sent: Thursday, June 22, 2006 4:59 PM To: libmesh-users@... Subject: [Libmesh-users] Question about Dirichlet BC imposed by penaltymethod Hi, I know there are already many examples at libmesh web site using penalty method to impose Dirichlet BC. However none of of them is generalized eigenvalue problem related and most of them have vector as the RHS. My problem is to solve the eigenvalue of a stationary 3D Schrodinger equation. The element loop for the mass and stiffness matrices is the following: ------------------------------------------------------------------------ ------------------------------------------------------------------ for (unsigned int qp=3D0; qpn_sides(); side++) if (elem->neighbor(side) =3D=3D NULL) { const std::vector >& phi_face =3D fe_face->get_phi(); const std::vector >& dphi_face =3D fe_face->get_dphi(); const std::vector& JxW_face =3D = fe_face->get_JxW(); const std::vector& qface_point =3D fe_face->get_xyz(); fe_face->reinit(elem, side); for (unsigned int qp=3D0; qp

 Re: [Libmesh-users] Question about Dirichlet BC imposed by penaltymethod From: Kirk, Benjamin \(JSC-EG\) - 2006-06-23 17:18:40 ```What you want is the linear system to reduce to the L2 projection of the boundary data for the boundary degrees of freedom. The way we tend to do that is to add the L2 projection, scaled by some *huge* penalty value, to the matrix and RHS. In floating point arithmetic the idea is that you add these huge values to the matrix and the other values in the matrix are no longer significant and essentially are 0. One reason for doing this is that you only need to assess boundary conditions in elements that have a full DIM-1 face on the boundary and can ignore elements with lower-dimensional boundary traces (edges in 3D, nodes in 2D/3D). =20 I am not familiar with your formulation, but I am almost certain you want the penalty value multiplying the Me() term and not the Ke() term. Again, the notion is to add the L2 projection of the boundary data, and the mass matrix contains the terms necessary for the L2 projection on the matrix side. -Ben =20 -----Original Message----- From: libmesh-users-bounces@... [mailto:libmesh-users-bounces@...] On Behalf Of David Xu Sent: Thursday, June 22, 2006 4:59 PM To: libmesh-users@... Subject: [Libmesh-users] Question about Dirichlet BC imposed by penaltymethod Hi, I know there are already many examples at libmesh web site using penalty method to impose Dirichlet BC. However none of of them is generalized eigenvalue problem related and most of them have vector as the RHS. My problem is to solve the eigenvalue of a stationary 3D Schrodinger equation. The element loop for the mass and stiffness matrices is the following: ------------------------------------------------------------------------ ------------------------------------------------------------------ for (unsigned int qp=3D0; qpn_sides(); side++) if (elem->neighbor(side) =3D=3D NULL) { const std::vector >& phi_face =3D fe_face->get_phi(); const std::vector >& dphi_face =3D fe_face->get_dphi(); const std::vector& JxW_face =3D = fe_face->get_JxW(); const std::vector& qface_point =3D fe_face->get_xyz(); fe_face->reinit(elem, side); for (unsigned int qp=3D0; qp
 Re: [Libmesh-users] Question about Dirichlet BC imposed by penaltymethod From: David Xu - 2006-06-23 22:53:56 ```Hi Ben, Thanks for your reply. I kinda understand the concept of the penalty scheme. The part I'm not sure is why we use L2 projection of the boundary data? Say, I just want to set the solution u=0 at the boudaries of a cubic domain and I have stiffness matrix Ke and mass matrix Me, the generalized eigenvalue problem equation is Ke*u=lamda*Me*u, the RHS is Me, not a vector. How would you impose penalty on Ke and Me to achieve u=0 BC? Some pseudo code would be even more helpful and highly appreciated. Regards, David On 6/23/06, Kirk, Benjamin (JSC-EG) wrote: > What you want is the linear system to reduce to the L2 projection of the > boundary data for the boundary degrees of freedom. The way we tend to > do that is to add the L2 projection, scaled by some *huge* penalty > value, to the matrix and RHS. In floating point arithmetic the idea is > that you add these huge values to the matrix and the other values in the > matrix are no longer significant and essentially are 0. > > One reason for doing this is that you only need to assess boundary > conditions in elements that have a full DIM-1 face on the boundary and > can ignore elements with lower-dimensional boundary traces (edges in 3D, > nodes in 2D/3D). > > I am not familiar with your formulation, but I am almost certain you > want the penalty value multiplying the Me() term and not the Ke() term. > Again, the notion is to add the L2 projection of the boundary data, and > the mass matrix contains the terms necessary for the L2 projection on > the matrix side. > > -Ben > > > > > -----Original Message----- > From: libmesh-users-bounces@... > [mailto:libmesh-users-bounces@...] On Behalf Of David > Xu > Sent: Thursday, June 22, 2006 4:59 PM > To: libmesh-users@... > Subject: [Libmesh-users] Question about Dirichlet BC imposed by > penaltymethod > > Hi, > > I know there are already many examples at libmesh web site using penalty > method to impose Dirichlet BC. However none of of them is generalized > eigenvalue problem related and most of them have vector as the RHS. My > problem is to solve the eigenvalue of a stationary 3D Schrodinger > equation. > > The element loop for the mass and stiffness matrices is the following: > > ------------------------------------------------------------------------ > ------------------------------------------------------------------ > for (unsigned int qp=0; qp { > const Real x = q_point[qp](0); > const Real y = q_point[qp](1); > const Real z = q_point[qp](2); > for (unsigned int i=0; i for (unsigned int j=0; j { > Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]; > Ke(i,j) += JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] + > (.5*pow(x, 2.)+.72*pow(y, 2.)+.845*pow(z, 2.))*phi[i][qp]*phi[j][qp]); > } > } > ------------------------------------------------------------------------ > ------------------------------------------------------------------ > > How do I impose the simple Dirichlet BC( u=0 at the boundary) using > penalty method? Does the following code look right? > > ------------------------------------------------------------------------ > ------------------------------------------------------------------ > for (unsigned int side=0; siden_sides(); side++) > if (elem->neighbor(side) == NULL) > { > const std::vector >& phi_face = > fe_face->get_phi(); > const std::vector >& > dphi_face = fe_face->get_dphi(); > const std::vector& JxW_face = fe_face->get_JxW(); > const std::vector& qface_point = > fe_face->get_xyz(); > fe_face->reinit(elem, side); > > for (unsigned int qp=0; qp { > const Real xf = qface_point[qp](0); > const Real yf = qface_point[qp](1); > const Real zf = qface_point[qp](2); > > const Real penalty = 1.e10; > > for (unsigned int i=0; i for (unsigned int j=0; > j { > Me(i,j) += > JxW_face[qp]*0*penalty*phi_face[i][qp]*phi_face[j][qp]; > Ke(i,j) += > JxW_face[qp]*penalty*(.5*dphi_face[i][qp]*dphi_face[j][qp] + (.5*pow(xf, > 2.)+.72*pow(yf, 2.)+.845*pow(zf, 2.))*phi_face[i][qp]*phi_face[j][qp]); > } > } > } > ------------------------------------------------------------------------ > ------------------------------------------------------------------ > > Thanks a lot in advance! > > David > > Using Tomcat but need to do more? Need to support web services, > security? > Get stuff done quickly with pre-integrated technology to make your job > easier Download IBM WebSphere Application Server v.1.0.1 based on Apache > Geronimo > http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642 > _______________________________________________ > Libmesh-users mailing list > Libmesh-users@... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > ```