You can subscribe to this list here.
2003 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}
(2) 
_{Oct}
(2) 
_{Nov}
(27) 
_{Dec}
(31) 

2004 
_{Jan}
(6) 
_{Feb}
(15) 
_{Mar}
(33) 
_{Apr}
(10) 
_{May}
(46) 
_{Jun}
(11) 
_{Jul}
(21) 
_{Aug}
(15) 
_{Sep}
(13) 
_{Oct}
(23) 
_{Nov}
(1) 
_{Dec}
(8) 
2005 
_{Jan}
(27) 
_{Feb}
(57) 
_{Mar}
(86) 
_{Apr}
(23) 
_{May}
(37) 
_{Jun}
(34) 
_{Jul}
(24) 
_{Aug}
(17) 
_{Sep}
(50) 
_{Oct}
(24) 
_{Nov}
(10) 
_{Dec}
(60) 
2006 
_{Jan}
(47) 
_{Feb}
(46) 
_{Mar}
(127) 
_{Apr}
(19) 
_{May}
(26) 
_{Jun}
(62) 
_{Jul}
(47) 
_{Aug}
(51) 
_{Sep}
(61) 
_{Oct}
(42) 
_{Nov}
(50) 
_{Dec}
(33) 
2007 
_{Jan}
(60) 
_{Feb}
(55) 
_{Mar}
(77) 
_{Apr}
(102) 
_{May}
(82) 
_{Jun}
(102) 
_{Jul}
(169) 
_{Aug}
(117) 
_{Sep}
(80) 
_{Oct}
(37) 
_{Nov}
(51) 
_{Dec}
(43) 
2008 
_{Jan}
(71) 
_{Feb}
(94) 
_{Mar}
(98) 
_{Apr}
(125) 
_{May}
(54) 
_{Jun}
(119) 
_{Jul}
(60) 
_{Aug}
(111) 
_{Sep}
(118) 
_{Oct}
(125) 
_{Nov}
(119) 
_{Dec}
(94) 
2009 
_{Jan}
(109) 
_{Feb}
(38) 
_{Mar}
(93) 
_{Apr}
(88) 
_{May}
(29) 
_{Jun}
(57) 
_{Jul}
(53) 
_{Aug}
(48) 
_{Sep}
(68) 
_{Oct}
(151) 
_{Nov}
(23) 
_{Dec}
(35) 
2010 
_{Jan}
(84) 
_{Feb}
(60) 
_{Mar}
(184) 
_{Apr}
(112) 
_{May}
(60) 
_{Jun}
(90) 
_{Jul}
(23) 
_{Aug}
(70) 
_{Sep}
(119) 
_{Oct}
(27) 
_{Nov}
(47) 
_{Dec}
(54) 
2011 
_{Jan}
(22) 
_{Feb}
(19) 
_{Mar}
(92) 
_{Apr}
(93) 
_{May}
(35) 
_{Jun}
(91) 
_{Jul}
(32) 
_{Aug}
(61) 
_{Sep}
(7) 
_{Oct}
(69) 
_{Nov}
(81) 
_{Dec}
(23) 
2012 
_{Jan}
(64) 
_{Feb}
(95) 
_{Mar}
(35) 
_{Apr}
(36) 
_{May}
(63) 
_{Jun}
(98) 
_{Jul}
(70) 
_{Aug}
(171) 
_{Sep}
(149) 
_{Oct}
(64) 
_{Nov}
(67) 
_{Dec}
(126) 
2013 
_{Jan}
(108) 
_{Feb}
(104) 
_{Mar}
(171) 
_{Apr}
(133) 
_{May}
(108) 
_{Jun}
(100) 
_{Jul}
(93) 
_{Aug}
(126) 
_{Sep}
(74) 
_{Oct}
(59) 
_{Nov}
(145) 
_{Dec}
(93) 
2014 
_{Jan}
(38) 
_{Feb}
(45) 
_{Mar}
(26) 
_{Apr}
(41) 
_{May}
(125) 
_{Jun}
(70) 
_{Jul}
(61) 
_{Aug}
(66) 
_{Sep}
(60) 
_{Oct}
(110) 
_{Nov}
(27) 
_{Dec}
(30) 
2015 
_{Jan}
(43) 
_{Feb}
(67) 
_{Mar}
(71) 
_{Apr}
(92) 
_{May}
(39) 
_{Jun}
(15) 
_{Jul}
(46) 
_{Aug}
(63) 
_{Sep}
(84) 
_{Oct}
(82) 
_{Nov}
(69) 
_{Dec}
(45) 
2016 
_{Jan}
(92) 
_{Feb}
(91) 
_{Mar}
(148) 
_{Apr}
(43) 
_{May}
(58) 
_{Jun}
(117) 
_{Jul}
(92) 
_{Aug}
(140) 
_{Sep}
(49) 
_{Oct}
(33) 
_{Nov}
(85) 
_{Dec}
(40) 
2017 
_{Jan}
(41) 
_{Feb}
(36) 
_{Mar}
(49) 
_{Apr}
(40) 
_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 







1

2

3

4
(1) 
5
(1) 
6

7
(8) 
8

9

10
(5) 
11
(6) 
12
(2) 
13
(1) 
14
(1) 
15
(1) 
16

17

18

19

20

21

22

23

24

25
(1) 
26

27

28

29

30

31






From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050125 21:30:11

I received an email a week or so ago, and for the life of me I can't find it... It asked specifically about Cygwin support. So, I'm sending this to the devel and users lists in hopes it finds the interested parties. I just installed cygwin on my windows machine (actually, cygwin is running under windows which is running in vmware under linux, talk about contrived!) here and was able to run the examples with no problem. Note, however, that shared libraries are not supported under cygwin at the current time. The following should work: ./configure disableshared make run_examples or at least it did for me with a fresh cygwin install using gcc3.3.3. Perhaps in the future I'll change configure so that it recogizes cygwin and disables shared libraries automatically.  Benjamin S. Kirk benjamin.kirk@... NASA/JSC, Applied Aerosciences & CFD Branch 2814839491 
From: Benjamin S. Kirk <benkirk@cf...>  20050115 03:54:44

I have used DG in libMesh for the Euler and NavierStokes equations in compressible gas dynamics. The MONOMIAL and XYZ finite element families are both discontinuous, in the first case the element shape functions are defined in terms of the local, element coordinates (xi,eta,zeta), whereas in the second case the shape functions are defined in terms of the global (x,y,z) coordinates. For DG it is especially useful to use the XYZ family. This makes the (many) face evaluations fairly cheap. As Roy points out, however, currently libMesh is uniform p. For the case of discontinuous elements it would be pretty straightforward to extend the code to variable p. I have just not done that yet because I figured we would tackle hp refinement for the continuous case and hash out all the details then. If hrefinement is insufficient for your case please let me know. Like I said, adding p refinement for discontinuous elements is pretty straightforward and would not require a lot of work. Ben Roy Stogner wrote: > On Wed, 12 Jan 2005, Aleksander Grm wrote: > >> I have some strange illposed hyperbolic PDE. I hope that I could >> obtain solutions with discontinuous galerkin FEM with hprefinement. >> I saw that libMesh has very powerful things inside and I was wondering >> if libMesh can help me with that. >> Can you help me with answer how much or what do I have to upgrade or >> how far can I do it with libMesh now. > > > I'm sure one of the primary developers would have a better answer, but > I'll take a crack at it: > > The h refinement is no problem; libMesh is excellent for hangingnode > adaptive h refinement. I'll bet DG complicates this with the element > jump terms, but those complications should be addressible in your > system assembly code without any internal library changes. > > Using DG methods isn't a big problem  I don't know if you can do it > with the library as is (since global continuity is enforced by the way > the individual finite element classes associate degrees of freedom > with edge and corner nodes), but you could probably create a working > discontinuous finite element by just copying the C0 hierarchic element > and telling your new version to associate every degree of freedom with > the element interior. > > The hardest part of the problem here seems to be the adaptive p > refinement. LibMesh currently assumes that each scalar variable in > your system is associated with a single finite element object and each > finite element object has the same number of degrees of freedom on > every cell. I'm not sure if you can get around that without either > leaving tons of inefficient unused DoFs in your global system or > digging deep into the libMesh guts. You could combine uniform p > refinement with adaptive h refinement pretty easily, but I'm sure you > want more than that. >  > Roy Stogner > > >  > 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: Roy Stogner <roystgnr@ic...>  20050114 02:07:17

On Wed, 12 Jan 2005, Aleksander Grm wrote: > I have some strange illposed hyperbolic PDE. I hope that I could obtain > solutions with discontinuous galerkin FEM with hprefinement. > I saw that libMesh has very powerful things inside and I was wondering if > libMesh can help me with that. > Can you help me with answer how much or what do I have to upgrade or how far > can I do it with libMesh now. I'm sure one of the primary developers would have a better answer, but I'll take a crack at it: The h refinement is no problem; libMesh is excellent for hangingnode adaptive h refinement. I'll bet DG complicates this with the element jump terms, but those complications should be addressible in your system assembly code without any internal library changes. Using DG methods isn't a big problem  I don't know if you can do it with the library as is (since global continuity is enforced by the way the individual finite element classes associate degrees of freedom with edge and corner nodes), but you could probably create a working discontinuous finite element by just copying the C0 hierarchic element and telling your new version to associate every degree of freedom with the element interior. The hardest part of the problem here seems to be the adaptive p refinement. LibMesh currently assumes that each scalar variable in your system is associated with a single finite element object and each finite element object has the same number of degrees of freedom on every cell. I'm not sure if you can get around that without either leaving tons of inefficient unused DoFs in your global system or digging deep into the libMesh guts. You could combine uniform p refinement with adaptive h refinement pretty easily, but I'm sure you want more than that.  Roy Stogner 
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: Aleksander Grm <agrm@us...>  20050112 14:15:18

Dear all, I have some strange illposed hyperbolic PDE. I hope that I could obtain solutions with discontinuous galerkin FEM with hprefinement. I saw that libMesh has very powerful things inside and I was wondering if libMesh can help me with that. Can you help me with answer how much or what do I have to upgrade or how far can I do it with libMesh now. I'm aware that I should do a lot of work, but libMesh seems the right answer. Many thanks for answer and help. Aleksander 
From: Aleksander Grm <grm@ma...>  20050112 11:42:04

Dear all, I have some strange illposed hyperbolic PDE. I hope that I could obtain solutions with discontinuous galerkin FEM with hprefinement. I saw that libMesh has very powerful things inside and I was wondering if libMesh can help me with that. Can you help me with answer how much or what do I have to upgrade or how far can I do it with libMesh now. I'm aware that I should do a lot of work, but libMesh seems the right answer. Many thanks for answer and help. Aleksander 
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...>  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 14:40:14

Martin L=FCthi writes: > Hi >=20 > Maybe a stupid question: Why is the following construct used in all > the examples to determine that an element lives on the boundary? >=20 > for (unsigned int s=3D0; s<elem>n_sides(); s++) > if (elem>neighbor(s) =3D=3D NULL) >=20 > Would'nt this be faster / more obvious if one would use >=20 > if (elem>on_boundary()) >=20 > and then use the boundary side iterators (Elem::boundary_sides_begin= )? Heh, yes probably. The boundary side iterators were only recently added, the rest of the code hasn't been updated yet. John 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050111 14:30:37

Absolutely. This is just a case of the examples predating this method = and the boundary side iterators. They should be updated. Along those lines, I have a new approach for building sides that is = less expensive in terms of memory allocation and data copying which will be included in libMesh0.5.0 realsoonnow. Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Martin = L=FCthi Sent: Monday, January 10, 2005 7:28 PM To: libmeshusers@... Subject: [Libmeshusers] Testing if element is on boundary Hi Maybe a stupid question: Why is the following construct used in all the examples to determine that an element lives on the boundary? for (unsigned int s=3D0; s<elem>n_sides(); s++) if (elem>neighbor(s) =3D=3D NULL) Would'nt this be faster / more obvious if one would use if (elem>on_boundary()) and then use the boundary side iterators (Elem::boundary_sides_begin)? Thanks for clarification Martin  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: <luthi@gi...>  20050111 01:22:32

Hi Maybe a stupid question: Why is the following construct used in all the examples to determine that an element lives on the boundary? for (unsigned int s=0; s<elem>n_sides(); s++) if (elem>neighbor(s) == NULL) Would'nt this be faster / more obvious if one would use if (elem>on_boundary()) and then use the boundary side iterators (Elem::boundary_sides_begin)? Thanks for clarification Martin 
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...>  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: 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: 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 
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: Steffen Petersen <steffen.petersen@tu...>  20050107 22:21:29

Ups, I just saw that a different connectivity is used for gmv output. Hence, element one, three, five, and six are likely to cause a negative Jacobian. Steffen Pollini Andrea wrote: > On Friday 07 January 2005 15:27, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > >>A negative Jacobian occurs when you have an overlydistorted or inverted >>element. It can also be caused by having elements numbered in a >>nonconsistent manner, e.g. clockwise instead of counterclockwise. >> >>Where did your mesh come from? If it is small, perhaps you can post it and >>I can take a look? >> >> > > > Hi Ben, > the mesh is a cube divided into 6 tetrahedra, and I have build it by > hand... after I uniformly refine it. > > the mesh is below... > > > > gmvinput ascii > > nodes 8 > 0 0 1 1 0 0 1 1 > 0 1 0 1 0 1 0 1 > 0 0 0 0 1 1 1 1 > > cells 6 > tet 4 > 1 4 8 3 > tet 4 > 1 7 8 3 > tet 4 > 1 7 8 5 > tet 4 > 1 4 8 2 > tet 4 > 1 6 8 2 > tet 4 > 1 6 8 5 > > material 1 0 > proc_0 > 1 > 1 > 1 > 1 > 1 > 1 > > > endgmv > > 
From: Steffen Petersen <steffen.petersen@tu...>  20050107 21:03:41

You can check your element definitions with the documentation of class Tet4. All tetrahedral elements have to be defined in this way: TET4: 3 o /\ /  \ /  \ 0 o......o 2 \  / \  / \/ o 1 In the example mesh you sent below it looks as if e.g. the second element will lead to a negative Jacobian. Steffen Pollini Andrea wrote: > On Friday 07 January 2005 15:27, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > >>A negative Jacobian occurs when you have an overlydistorted or inverted >>element. It can also be caused by having elements numbered in a >>nonconsistent manner, e.g. clockwise instead of counterclockwise. >> >>Where did your mesh come from? If it is small, perhaps you can post it and >>I can take a look? >> >> > > > Hi Ben, > the mesh is a cube divided into 6 tetrahedra, and I have build it by > hand... after I uniformly refine it. > > the mesh is below... > > > > gmvinput ascii > > nodes 8 > 0 0 1 1 0 0 1 1 > 0 1 0 1 0 1 0 1 > 0 0 0 0 1 1 1 1 > > cells 6 > tet 4 > 1 4 8 3 > tet 4 > 1 7 8 3 > tet 4 > 1 7 8 5 > tet 4 > 1 4 8 2 > tet 4 > 1 6 8 2 > tet 4 > 1 6 8 5 > > material 1 0 > proc_0 > 1 > 1 > 1 > 1 > 1 > 1 > > > endgmv > > 
From: Pollini Andrea <pollini@dm...>  20050107 17:06:37

On Friday 07 January 2005 15:27, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > A negative Jacobian occurs when you have an overlydistorted or inverted > element. It can also be caused by having elements numbered in a > nonconsistent manner, e.g. clockwise instead of counterclockwise. > > Where did your mesh come from? If it is small, perhaps you can post it and > I can take a look? > > Hi Ben, the mesh is a cube divided into 6 tetrahedra, and I have build it by hand... after I uniformly refine it. the mesh is below... gmvinput ascii nodes 8 0 0 1 1 0 0 1 1 0 1 0 1 0 1 0 1 0 0 0 0 1 1 1 1 cells 6 tet 4 1 4 8 3 tet 4 1 7 8 3 tet 4 1 7 8 5 tet 4 1 4 8 2 tet 4 1 6 8 2 tet 4 1 6 8 5 material 1 0 proc_0 1 1 1 1 1 1 endgmv thanks in advance, Pollini Andrea 
From: Michael Schindler <mschindler@us...>  20050107 15:05:06

Hello Ben, On 07.01.05, KIRK, BENJAMIN (JSCEG) (NASA) wrote: > For example, consider Laplace's problem in 2D: > >  div(grad(u)) = 0 > > with u=g on some part of the boundary, du/dn=h on the remainder. > > In the weak statement you perform integrationbyparts, and the term du/dn > appears in a boundary integral and you can use this to impose a Neumann BC > in the RHS of the system. See, for example, > http://cfdlab.ae.utexas.edu/~benkirk/seminar/talk.pdf, Especially page 4. Do you mean just to replace du/dn in the boundary integral by h in order to "enforce" du/dn = h ?? This would be the FEM standard procedure, I guess. Exactly at this point I have some doubts. Consider a different setup, where I add the equation du/dn = h weighted with some penalty, to the linear system to be solved. This works quite well. The first prodecure is equivalent to this penaltyprocedure, if one sets the penalty to 1 which is not a really good value for a penalty. > It's not clear that you need to use the constraint matrix for this case... Therefore, I wanted to look for an alternative approach to enforce the boundary condition  and came across the constraint matrices. > For implementation, you simply need to define a finite element object that > lives on the boundary and can be used to integrate along the edge to provide > the required term. Clear? Yes, clear. For an element I would get a constraining equation (again for the case dn/du = h) h = \sum_i w_i u_i where u_i are the degrees of freedom and w_i are some weights determined by the geometry of the element. Now comes the next difficulty. I have an idea how the method DofMap::constrain_matrix_and_vector works. It creates a matrix C from the constraints and uses these constraints while inserting the correct values into the big system matrix. I would expect this _only_ to work for h = 0. Do I miss something here? > Unfortunately, none of the examples show this. If you have any more > questions let me know, it would be straightforward to modify ex14 to show > this procedure since it has the exact solution's derivative available. > Regardless, we should probably do that since it seems silly to have 14 > examples, all with Dirichlet BCs! This would be great! Thanks, Michael.  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 
From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050107 14:43:34

It's not clear that you need to use the constraint matrix for this = case... For example, consider Laplace's problem in 2D: =20  div(grad(u)) =3D 0 with u=3Dg on some part of the boundary, du/dn=3Dh on the remainder. In the weak statement you perform integrationbyparts, and the term = du/dn appears in a boundary integral and you can use this to impose a Neumann = BC in the RHS of the system. See, for example, http://cfdlab.ae.utexas.edu/~benkirk/seminar/talk.pdf, Especially page = 4. In a similar way, Robin BCs can be implemented by letting h =3D alpha u = + beta du/dn. In this case you have two boundary terms that contribute to the system matrix. For implementation, you simply need to define a finite element object = that lives on the boundary and can be used to integrate along the edge to = provide the required term. Clear? Unfortunately, none of the examples show this. If you have any more questions let me know, it would be straightforward to modify ex14 to = show this procedure since it has the exact solution's derivative available. Regardless, we should probably do that since it seems silly to have 14 examples, all with Dirichlet BCs! Ben =20 Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Michael Schindler Sent: Friday, January 07, 2005 3:26 AM To: libmeshusers@... Subject: [Libmeshusers] Constraint matrix for boundary conditions Hello, I would like to enforce nonDirichlet boundary conditions via the = constraint matrix. Does someone have an example how to implement this for e.g. = Neumann or Robin BC? Thanks, Michael. =20 "A mathematician is a device for turning coffee into theorems" Paul Erd=F6s.  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: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050107 14:27:05

A negative Jacobian occurs when you have an overlydistorted or inverted element. It can also be caused by having elements numbered in a nonconsistent manner, e.g. clockwise instead of counterclockwise. Where did your mesh come from? If it is small, perhaps you can post it and I can take a look? Ben Original Message From: libmeshusersadmin@... [mailto:libmeshusersadmin@...] On Behalf Of Pollini Andrea Sent: Friday, January 07, 2005 6:55 AM To: libmeshusers@... Subject: [Libmeshusers] problem with FEBase class Hi all, why when I call the method reinit in class FEBase I receive the error ERROR: negative Jacobian: 0.000244141 [0] src/fe/fe_map.C, line 405, compiled Jan 5 2005 at 11:40:52 /bin/sh: line 1: 11242 Abortito ./libmeshsteffan Press Enter to continue! thank in advance for every advice, Pollini Andrea 
From: Pollini Andrea <nelson@nu...>  20050107 12:56:48

Hi all, why when I call the method reinit in class FEBase I receive the error ERROR: negative Jacobian: 0.000244141 [0] src/fe/fe_map.C, line 405, compiled Jan 5 2005 at 11:40:52 /bin/sh: line 1: 11242 Abortito ./libmeshsteffan Press Enter to continue! thank in advance for every advice, Pollini Andrea 
From: Michael Schindler <mschindler@us...>  20050107 09:25:42

Hello, I would like to enforce nonDirichlet boundary conditions via the constraint matrix. Does someone have an example how to implement this for e.g. Neumann or Robin BC? Thanks, Michael.  "A mathematician is a device for turning coffee into theorems" Paul Erdös. 