From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20060623 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 DIM1 face on the boundary and can ignore elements with lowerdimensional 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: 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=3D0; qp<qrule.n_points(); qp++) { const Real x =3D q_point[qp](0); const Real y =3D q_point[qp](1); const Real z =3D q_point[qp](2); for (unsigned int i=3D0; i<phi.size(); i++) for (unsigned int j=3D0; j<phi.size(); j++) { Me(i,j) +=3D JxW[qp]*phi[i][qp]*phi[j][qp]; Ke(i,j) +=3D 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=3D0 at the boundary) using penalty method? Does the following code look right?   for (unsigned int side=3D0; side<elem>n_sides(); side++) if (elem>neighbor(side) =3D=3D NULL) { const std::vector<std::vector<Real> >& phi_face =3D fe_face>get_phi(); const std::vector<std::vector<RealGradient> >& dphi_face =3D fe_face>get_dphi(); const std::vector<Real>& JxW_face =3D = fe_face>get_JxW(); const std::vector<Point >& qface_point =3D fe_face>get_xyz(); fe_face>reinit(elem, side); for (unsigned int qp=3D0; qp<qface.n_points(); qp++) { const Real xf =3D qface_point[qp](0); const Real yf =3D qface_point[qp](1); const Real zf =3D qface_point[qp](2); const Real penalty =3D 1.e10; for (unsigned int i=3D0; i<phi_face.size(); = i++) for (unsigned int j=3D0; j<phi_face.size(); j++) { Me(i,j) +=3D JxW_face[qp]*0*penalty*phi_face[i][qp]*phi_face[j][qp]; Ke(i,j) +=3D 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=3Dlnk&kid=3D120709&bid=3D263057&dat=3D= 121642 _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 