From: John Peterson <peterson@cf...>  20050111 00:30:07

Hi Mark, Can you post a compilable version of your code? My only other guess is you could have a sign error since it's going to the right hand side. J Mark Blome writes: > > Hi Michael, > > I am using libmesh to solve for the electric potential of a current source > within a 3d cubic domain. The calculations are used for a geoelectrical > forward solver. On the top side of the domain, I apply a Neuman boundary > condition (du/dn=0) and on the other sides I use Dirichlet boundary > conditions (u=0). > The second term > div( (sigmasigma0) grad Un) > is added to remove the singularity arising due to the point current source. > When I discretize the second term with > > Fe(i) +=  (sigmasigma0) * JxW[qp] * (gradVn[qp] * dphi[i][qp]) > (of course, as John noticed, I forgot to mention the surface integral term, > but in my case it is zero anyway cause sigma=sigma0 at the outer boundary. ) > > I get wrong results: The resulting potential U is the potential that does not > contain the singular potential Un=1/(2*pi) *1/r anymore. > For the case of a homogenous cube with a high conducting cube centered inside, > the introduction of the discretization above results in a strange dipole > field within the high conducting cube, which physically doesnt make any > sense. > Thanks for the help so far, > > Mark > > > > Am Montag 10 Januar 2005 19:28 schrieb Michael Schindler: > > Hi Mark, > > > > On 10.01.05, Mark Blome wrote: > > > Hi everybody, > > > > > > I am trying to descretize a laplace equation which contains a second > > > term. In my opinion the problem should be very easy, but I become more > > > and more confused the more I think about it. The laplace equation reads: > > > div(sigma grad U) + div( (sigmasigma0) grad Un) =0 > > > where Un is a potential for which the equation is known ( Un=j0/(2 * PI * > > > sigma0) * 1/r), hence the term has to go to the right hand side. > > > > > > For the term div(sigma grad U)=0 I simply get matrix terms: > > > Ke(i,j) +=  sigma * JxW[qp]*(dphi[i][qp]*dphi[j][qp]); > > > But what to add for div( (sigmasigma0) grad Un)? > > > > The FEM method is to find a weak solution with a testfunction phi > > > > int div(sigma grad U) phi = int div( (sigmasigma0) grad Un) phi > > > > Depending on sigma you can use the Galerkinansatz where phi is the > > same function you approximate U with. If sigma is constant, omitting > > the surface terms. > > > >  sigma int grad(U)*grad(phi) =  (sigmasigma0) int grad(Un)*grad(phi) > > > > This is what you tried first. > > > > > I tried: > > > Fe(i) +=  (sigmasigma0) * JxW[qp] * (gradVn[qp] * dphi[i][qp]) > > > ; where gradVn[qp] simply is grad(Un) analytically calculated at the > > > quadrature points. But it doesnt work like this. > > > > What actually does not work? What boundary conditions do you use? > > > > Michael. > > >  > The SF.Net email is sponsored by: Beat the postholiday blues > Get a FREE limited edition SourceForge.net tshirt from ThinkGeek. > It's fun and FREE  well, almost....http://www.thinkgeek.com/sfshirt > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers  John W. Peterson The University of Texas at Austin Web: http://www.cfdlab.ae.utexas.edu Phone: 5126196534 (H)  5124714069 (W) email: peterson@... 