From: Mark Blome <mblome@au...>  20050111 16:06:57

Hi John and Michael, thanks for your comments on this problem! I think I finally found the mistake. As usually in a place where I never expected it: The right hand side terms where correct, but for the calculation of the gradient of the potential Un I used something like double x= qPoints[qp](0)source.pos(0) ; double y= qPoints[qp](1)source.pos(1) ; double z= qPoints[qp](2)source.pos(2) ; double r=sqrt( pow(x,2) + pow(y,2) + pow(z,2) ); double cpart=  source.j0/(2*libMesh::pi * sigma0); RealGradient gradVn; gradVn(0)=cpart*x/pow(r,3); gradVn(1)=cpart*y/pow(r,3); gradVn(2)=cpart*z/pow(r,3); in the quadrature loop. (source.pos is the coordinate of the current source and the qpoints are from the quad rule: const std::vector<Point> qPoints=qrule.get_points();) I was simply not aware that "get_points()" returns local coordinates on a reference element and not global coodinates. If I use eg "elem>point(0)" instead of qPoints[qp] the results are fine! Is there an easy way to convert the local reference element coordinates to global coordinates? Mark Am Dienstag 11 Januar 2005 01:29 schrieb John Peterson: > 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 