From: Benjamin S. K. <be...@cf...> - 2005-01-13 02:24:41
|
Ah, a simple mistake! You are right, the qrule.get_points() will return a reference to the quadrature points in local (element) coordinates. These are then mapped to physical space when you call fe.reinit(elem) or fe.reinit(elem,side). The quadrature points in physical space can be obtained from the FE object: const std::vector<Point>& q_point = fe->get_xyz(); See, for example, line 266 in ex4.C. -Ben Mark Blome wrote: > 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( (sigma-sigma0) grad Un) >> > is added to remove the singularity arising due to the point current >> > source. When I discretize the second term with >> > >> > Fe(i) += - (sigma-sigma0) * 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( (sigma-sigma0) 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( (sigma-sigma0) grad Un)? >> > > >> > > The FEM method is to find a weak solution with a test-function phi >> > > >> > > int div(sigma grad U) phi = int div( (sigma-sigma0) grad Un) phi >> > > >> > > Depending on sigma you can use the Galerkin-ansatz where phi is the >> > > same function you approximate U with. If sigma is constant, omitting >> > > the surface terms. >> > > >> > > - sigma int grad(U)*grad(phi) = - (sigma-sigma0) int >> > > grad(Un)*grad(phi) >> > > >> > > This is what you tried first. >> > > >> > > > I tried: >> > > > Fe(i) += - (sigma-sigma0) * 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 post-holiday blues >> > Get a FREE limited edition SourceForge.net t-shirt from ThinkGeek. >> > It's fun and FREE -- well, almost....http://www.thinkgeek.com/sfshirt >> > _______________________________________________ >> > Libmesh-users mailing list >> > Lib...@li... >> > https://lists.sourceforge.net/lists/listinfo/libmesh-users > > > > ------------------------------------------------------- > The SF.Net email is sponsored by: Beat the post-holiday blues > Get a FREE limited edition SourceForge.net t-shirt from ThinkGeek. > It's fun and FREE -- well, almost....http://www.thinkgeek.com/sfshirt > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users |