From: David Xu <dxu@my...>  20060623 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 (JSCEG) <benjamin.kirk1@...> 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 DIM1 face on the boundary and > can ignore elements with lowerdimensional 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: libmeshusersbounces@... > [mailto:libmeshusersbounces@...] On Behalf Of David > Xu > Sent: Thursday, June 22, 2006 4:59 PM > To: libmeshusers@... > Subject: [Libmeshusers] 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<qrule.n_points(); 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<phi.size(); i++) > for (unsigned int j=0; j<phi.size(); 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; side<elem>n_sides(); side++) > if (elem>neighbor(side) == NULL) > { > const std::vector<std::vector<Real> >& phi_face = > fe_face>get_phi(); > const std::vector<std::vector<RealGradient> >& > dphi_face = fe_face>get_dphi(); > const std::vector<Real>& JxW_face = fe_face>get_JxW(); > const std::vector<Point >& qface_point = > fe_face>get_xyz(); > fe_face>reinit(elem, side); > > for (unsigned int qp=0; qp<qface.n_points(); 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<phi_face.size(); i++) > for (unsigned int j=0; > j<phi_face.size(); 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 preintegrated technology to make your job > easier Download IBM WebSphere Application Server v.1.0.1 based on Apache > Geronimo > http://sel.asus.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642 > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 