From: John Peterson <peterson@cf...>  20050111 16:16:29

Mark Blome writes: > > 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? Yes, use get_xyz(). John 
From: Mark Blome <mblome@au...>  20050110 11:05:32

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)? 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. Then I tried to express Un using the form functions, ie Un=sum ( Un_0 * phi_i ), which would result to Fe(i) =  (sigmasigma0) * JxW[qp] * ( (dphi[i][qp] * Vn[qp] + gradVn[qp] *phi[i][qp]) * dphi[i][qp] ) That doesnt work as well. I also tried Fe(i) = (sigmasigma0) * JxW[qp] * (dphi[i][qp] * dphi[i][qp]) * Vn[i]; which doesnt make so much sense to me and as I expected doesnt work eather. I think this is really an easy thing, but my thoughts seem to go in a completely wrong direction. I would be very happy for any comments on this, Mark 
From: Michael Schindler <mschindler@us...>  20050110 18:28:34

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.  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: Mark Blome <mblome@au...>  20050110 21:10:25

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. 
From: Mark Blome <mblome@au...>  20050110 21:10:25

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. 
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@... 
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 
From: John Peterson <peterson@cf...>  20050111 16:16:29

Mark Blome writes: > > 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? Yes, use get_xyz(). John 
From: Benjamin S. Kirk <benkirk@cf...>  20050113 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( (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 > > > >  > 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 
From: John Peterson <peterson@cf...>  20050110 15:02:26

Mark Blome writes: > > 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)? > 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. Hi. This is pretty close! The second term in your equation above, div( (sigmasigma0) grad Un) generates two terms when you multiply by a test function: \int_{Omega} div ( (sigmasigma0) grad(Un) ) * phi = + \int_{dOmega} (sigmasigma0) ( grad(Un) * n ) * phi dS (Term 1)  \int_{Omega} (sigmasigma0) grad(Un) * grad(phi) dx (Term 2) So you will have to code something pretty much like the first part, with the associated boundary term, but just put it on the rhs. John 