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}
(41) 
_{May}
(73) 
_{Jun}
(51) 
_{Jul}
(12) 
_{Aug}
(69) 
_{Sep}
(26) 
_{Oct}
(43) 
_{Nov}
(75) 
_{Dec}
(9) 
S  M  T  W  T  F  S 





1
(1) 
2
(2) 
3
(1) 
4

5
(1) 
6
(4) 
7
(2) 
8

9
(4) 
10
(1) 
11

12
(4) 
13
(4) 
14
(3) 
15
(4) 
16
(1) 
17

18
(2) 
19
(5) 
20
(1) 
21
(2) 
22
(1) 
23
(2) 
24
(1) 
25

26
(2) 
27
(2) 
28
(1) 
29
(6) 
30
(5) 

From: David Xu <dxu@my...>  20060630 19:09:50

On 6/29/06, John Peterson <peterson@...> wrote: > > Roy Stogner writes: > > On Thu, 29 Jun 2006, David Xu wrote: > > > > > When I was trying to generate 3D tetrahedral mesh using: > > > > > > MeshTools::Generation::build_cube (mesh, 20, 20, 20, 3., 3., 3., 3., > 3., > > > 3., TET10); > > > > > > I got: > > > > > > ERROR: Unrecognized 3D element type. > > > [0] src/mesh/mesh_generation.C, line 671, compiled Jun 13 2006 at > 23:07:35 > > > > > > I thought TET4 and TET10 elements were implemented in libmesh because > they > > > are listed under Generic 3D Finite Elements of libmesh Class DOCs > > > http://libmesh.sourceforge.net/doxygen/index.php > > > > Unfortunately, "implemented in libMesh" doesn't always mean > > "implemented in every single libMesh function". libMesh currently > > supports finite element calculations on tetrahedra with linear or > > quadratic Lagrange elements, but you need to be able to supply your > > own tet mesh (or get one out of tetgen, now). > > > > "Unrecognized 3D element type" is a misleading error message there; > > perhaps we should change it to something like "3D element type X is > > not supported by function Y". > > > > John Peterson has recently written some code that let him build cubes > > of tetrahedra, but I don't know if it's debugged and ready to add to > > CVS yet. > > I think TET4 and TET10 are supported by build_cube in the CVS version. > Note also that it works by creating 24 tets from each hex in an initial > mesh of (in your case) 20x20x20 hexes. That's probably way too many tets! > > my libmesh was built on the cvs version and I still got ERROR: Unrecognized 3D element type. I downloaded libmesh thru cvs on 6/12, was there any update since then? I'd like to try creating 24 tets from each hex, could you give me some sample code? Thanks, David 
From: Roy Stogner <roystgnr@ic...>  20060630 19:08:43

On Fri, 30 Jun 2006, David Knezevic wrote: > I've precomputed some (SparseMatrix<Number>) matrices that I can > reuse repeatedly in my computation, and I'd also like to be able to > compute the inverse of these matrices. I couldn't see any libMesh > functionality for computing matrix inverses, but I just wanted to > check I didn't overlook something. No, you didn't miss it, it's not there. The inverse of a sparse matrix is usually a dense matrix, and so you rarely want to compute it explicitly unless you're doing far more linear solves than you have matrix rows. You can trick PETSc into doing the equivalent of a matrix inverse by using pc_type lu, but you'll have to dig into the PETSc interface code to figure out how to get it to save that preconditioner between solves. > If there isn't anything in libMesh, would there be an easy way to do > this (perhaps using PETSc functionality)? > > Also, if I can compute these inverses, I'd then need to do a > SparseMatrix multiplication with a NumericVector. Again, I didn't see > multiplication defined anywhere, just wondering if multiplication is > already available or if I'd have to implement it myself? You're planning to multiply a matrix inverse by a vector  are you sure you don't just want to do a linear solve using that matrix and vector? You'll probably get better performance with a few ILU steps (saved between solves) and a few Krylov steps per solve, instead of one big Gaussian elimination. There is a NumericVector::add_vector() method that does U += A*V with a SparseMatrix A, though.  Roy 
From: David Knezevic <david.knezevic@ba...>  20060630 18:51:40

Hi all, I've precomputed some (SparseMatrix<Number>) matrices that I can reuse repeatedly in my computation, and I'd also like to be able to compute the inverse of these matrices. I couldn't see any libMesh functionality for computing matrix inverses, but I just wanted to check I didn't overlook something. If there isn't anything in libMesh, would there be an easy way to do this (perhaps using PETSc functionality)? Also, if I can compute these inverses, I'd then need to do a SparseMatrix multiplication with a NumericVector. Again, I didn't see multiplication defined anywhere, just wondering if multiplication is already available or if I'd have to implement it myself? Thanks, David 
From: li pan <li76pan@ya...>  20060630 07:15:12

Hi David, You could possiblly do like this. Use Tetgen to generate a cube of 6 tetelement. Save it in .node and .ele files. Load the files in you libmesh code and refine the cube with mesh_refinement.refine_uniformly(). see u pan Hi All, When I was trying to generate 3D tetrahedral mesh using: MeshTools::Generation::build_cube (mesh, 20, 20, 20, 3., 3., 3., 3., 3., 3., TET10); I got: ERROR: Unrecognized 3D element type. [0] src/mesh/mesh_generation.C, line 671, compiled Jun 13 2006 at 23:07:35 I thought TET4 and TET10 elements were implemented in libmesh because they are listed under Generic 3D Finite Elements of libmesh Class DOCs http://libmesh.sourceforge.net/doxygen/index.php Alternatively I tried tetgen with libmesh using: TetGenMeshInterface tetgen (mesh); tetgen.pointset_convexhull (); tetgen.triangulate_conformingDelaunayMesh(0., .005); And the program threw an error at equation_systems.init(); and [0] src/fe/fe_lagrange.C, line 614 Can anyone tell me what could be the problem? Thanks, David __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com 
From: John Peterson <peterson@cf...>  20060630 04:12:44

Roy Stogner writes: > On Thu, 29 Jun 2006, David Xu wrote: > > > When I was trying to generate 3D tetrahedral mesh using: > > > > MeshTools::Generation::build_cube (mesh, 20, 20, 20, 3., 3., 3., 3., 3., > > 3., TET10); > > > > I got: > > > > ERROR: Unrecognized 3D element type. > > [0] src/mesh/mesh_generation.C, line 671, compiled Jun 13 2006 at 23:07:35 > > > > I thought TET4 and TET10 elements were implemented in libmesh because they > > are listed under Generic 3D Finite Elements of libmesh Class DOCs > > http://libmesh.sourceforge.net/doxygen/index.php > > Unfortunately, "implemented in libMesh" doesn't always mean > "implemented in every single libMesh function". libMesh currently > supports finite element calculations on tetrahedra with linear or > quadratic Lagrange elements, but you need to be able to supply your > own tet mesh (or get one out of tetgen, now). > > "Unrecognized 3D element type" is a misleading error message there; > perhaps we should change it to something like "3D element type X is > not supported by function Y". > > John Peterson has recently written some code that let him build cubes > of tetrahedra, but I don't know if it's debugged and ready to add to > CVS yet. I think TET4 and TET10 are supported by build_cube in the CVS version. Note also that it works by creating 24 tets from each hex in an initial mesh of (in your case) 20x20x20 hexes. That's probably way too many tets! 
From: Roy Stogner <roystgnr@ic...>  20060629 22:08:04

On Thu, 29 Jun 2006, David Xu wrote: > When I was trying to generate 3D tetrahedral mesh using: > > MeshTools::Generation::build_cube (mesh, 20, 20, 20, 3., 3., 3., 3., 3., > 3., TET10); > > I got: > > ERROR: Unrecognized 3D element type. > [0] src/mesh/mesh_generation.C, line 671, compiled Jun 13 2006 at 23:07:35 > > I thought TET4 and TET10 elements were implemented in libmesh because they > are listed under Generic 3D Finite Elements of libmesh Class DOCs > http://libmesh.sourceforge.net/doxygen/index.php Unfortunately, "implemented in libMesh" doesn't always mean "implemented in every single libMesh function". libMesh currently supports finite element calculations on tetrahedra with linear or quadratic Lagrange elements, but you need to be able to supply your own tet mesh (or get one out of tetgen, now). "Unrecognized 3D element type" is a misleading error message there; perhaps we should change it to something like "3D element type X is not supported by function Y". John Peterson has recently written some code that let him build cubes of tetrahedra, but I don't know if it's debugged and ready to add to CVS yet. > Alternatively I tried tetgen with libmesh using: > > TetGenMeshInterface tetgen (mesh); > tetgen.pointset_convexhull (); > tetgen.triangulate_conformingDelaunayMesh(0., .005); > > And the program threw an error at > > equation_systems.init(); > > and [0] src/fe/fe_lagrange.C, line 614 > > Can anyone tell me what could be the problem? Well, I don't use tetgen so hopefully someone else can help you with that. If you're using the CVS libMesh, line 614 of fe_lagrange.C is the "you tried to use a polynomial order we can't handle" error() call. You'll have to get a full stack trace out of your debugger to find out where the offending code is, but some function in your program tried to use a quartic or higher Lagrange element, and on most geometric elements our Lagrange FE only handles linear and quadratic bases.  Roy Stogner 
From: Dong Xu <dongxu@sc...>  20060629 21:36:26

Hi All, When I was trying to generate 3D tetrahedral mesh using: MeshTools::Generation::build_cube (mesh, 20, 20, 20, 3., 3., 3., 3., 3., 3., TET10); I got: ERROR: Unrecognized 3D element type. [0] src/mesh/mesh_generation.C, line 671, compiled Jun 13 2006 at 23:07:35 I thought TET4 and TET10 elements were implemented in libmesh because they are listed under Generic 3D Finite Elements of libmesh Class DOCs http://libmesh.sourceforge.net/doxygen/index.php Alternatively I tried tetgen with libmesh using: TetGenMeshInterface tetgen (mesh); tetgen.pointset_convexhull(); tetgen.triangulate_conformingDelaunayMesh(0., .005); And the program threw an error at equation_systems.init(); and [0] src/fe/fe_lagrange.C, line 614 Can anyone tell me what could be the problem? Thanks, David 
From: David Xu <dxu@my...>  20060629 21:35:51

Hi All, When I was trying to generate 3D tetrahedral mesh using: MeshTools::Generation::build_cube (mesh, 20, 20, 20, 3., 3., 3., 3., 3., 3., TET10); I got: ERROR: Unrecognized 3D element type. [0] src/mesh/mesh_generation.C, line 671, compiled Jun 13 2006 at 23:07:35 I thought TET4 and TET10 elements were implemented in libmesh because they are listed under Generic 3D Finite Elements of libmesh Class DOCs http://libmesh.sourceforge.net/doxygen/index.php Alternatively I tried tetgen with libmesh using: TetGenMeshInterface tetgen (mesh); tetgen.pointset_convexhull (); tetgen.triangulate_conformingDelaunayMesh(0., .005); And the program threw an error at equation_systems.init(); and [0] src/fe/fe_lagrange.C, line 614 Can anyone tell me what could be the problem? Thanks, David 
From: David Xu <dxu@my...>  20060629 18:47:01

Thanks, John. I found that info in tetgen documentation too. I ran into another problem though. After I incorporated the tetgen sample code from Steffen, my program compiled ok. But at runtime, it stopped at the line: equation_systems.init(); and gave me the following error: [0] src/fe/fe_lagrange.C, line 614, compiled Jun 13 2006 at 22:52:46 Aborted I was wondering if the tetgen conversion from HEX8 to tetrahedral mesh generated some ElemType that is not recognized by libmesh? How did you guys use tetgen in your code? Thanks, David On 6/28/06, John Peterson <peterson@...> wrote: > > Quality constraint and area constraint. Take a look at the > documentation for tetgen if you want to know more. > > John > > David Xu writes: > > Hi All, > > > > What are the definitions for X, Y values in > > tetgen.triangulate_conformingDelaunayMesh(X, Y); > > > > Thanks, > > > > David > > > > > > On 6/21/06, Steffen Petersen <steffen.petersen@...> wrote: > > > > > > > > > Find attached an example file. > > > > > > Steffen > > > > > > > > > David Xu schrieb: > > > > Hi Steffen and John, > > > > > > > > I've been trying to figure out to how to generate tetrahedral mesh > in > > > > a simple cubic domain and I just happened to see your emails. I > don't > > > > know much about tetgen. Could you guys send me your files so that I > > > > can learn from them? > > > > > > > > Thanks a lot! > > > > > > > > David > > > > > > > > On 6/8/06, Steffen Petersen <steffen.petersen@...> wrote: > > > > > > > >> > > > >> I cannot see anything suspicious in the code. > > > >> I have just tested a simliar example which worked for me. > > > >> I will send the file to you. > > > >> > > > >> However, I have noticed some other problems. > > > >> In optimized mode I had to change the CXXFLAGS > > > >> in contrib/tetgen/Makefile > > > >> > > > >> #CXXFLAGS += DTRILIBRARY > > > >> CXXFLAGS := DNDEBUG DTRILIBRARY fPIC > > > >> > > > >> in order to get the example working. > > > >> > > > >> Steffen > > > >> > > > >> > > > >> > > > >> > Hi, > > > >> > > > > >> > I'm trying to use tetgen to generate a tetrahedral > > > >> > mesh in the cube [1,1]^3. > > > >> > > > > >> > So far, my code is the following: > > > >> > > > > >> > > > > >> > Mesh mesh(3); > > > >> > > > > >> > // add nodes at the eight vertices of the cube > > > >> > ... > > > >> > > > > >> > // Instantiate TetGen interface > > > >> > TetGenMeshInterface tetgen (mesh); > > > >> > > > > >> > > > > >> > // Determine the TRI3 elements which define the convex hull > > > >> > tetgen.pointset_convexhull(); > > > >> > > > > >> > // find neighbors, etc > > > >> > mesh.prepare_for_use(); > > > >> > > > > >> > > > > >> > // Generate tetrahedra > > > >> > tetgen.triangulate_conformingDelaunayMesh(X, Y); > > > >> > > > > >> > > > > >> > > > > >> > Now, here's where my trouble comes in: for any value of X and Y, > > > >> > I get the following error: > > > >> > > > > >> > tetgen.C:9469: long int tetgenmesh::flip(tetgenmesh::queue*, > > > >> > tetgenmesh::badface**): Assertion `fc != NONCONVEX' failed. > > > >> > > > > >> > > > > >> > What am I doing wrong? > > > >> > > > > >> > Thanks, > > > >> > John > > > >> > > > > >> > > 
From: Leonora Nixon <ohbhgj@fd...>  20060629 17:11:27

From: John Peterson <peterson@cf...>  20060629 03:30:55

Quality constraint and area constraint. Take a look at the documentation for tetgen if you want to know more. John David Xu writes: > Hi All, > > What are the definitions for X, Y values in > tetgen.triangulate_conformingDelaunayMesh(X, Y); > > Thanks, > > David > > > On 6/21/06, Steffen Petersen <steffen.petersen@...> wrote: > > > > > > Find attached an example file. > > > > Steffen > > > > > > David Xu schrieb: > > > Hi Steffen and John, > > > > > > I've been trying to figure out to how to generate tetrahedral mesh in > > > a simple cubic domain and I just happened to see your emails. I don't > > > know much about tetgen. Could you guys send me your files so that I > > > can learn from them? > > > > > > Thanks a lot! > > > > > > David > > > > > > On 6/8/06, Steffen Petersen <steffen.petersen@...> wrote: > > > > > >> > > >> I cannot see anything suspicious in the code. > > >> I have just tested a simliar example which worked for me. > > >> I will send the file to you. > > >> > > >> However, I have noticed some other problems. > > >> In optimized mode I had to change the CXXFLAGS > > >> in contrib/tetgen/Makefile > > >> > > >> #CXXFLAGS += DTRILIBRARY > > >> CXXFLAGS := DNDEBUG DTRILIBRARY fPIC > > >> > > >> in order to get the example working. > > >> > > >> Steffen > > >> > > >> > > >> > > >> > Hi, > > >> > > > >> > I'm trying to use tetgen to generate a tetrahedral > > >> > mesh in the cube [1,1]^3. > > >> > > > >> > So far, my code is the following: > > >> > > > >> > > > >> > Mesh mesh(3); > > >> > > > >> > // add nodes at the eight vertices of the cube > > >> > ... > > >> > > > >> > // Instantiate TetGen interface > > >> > TetGenMeshInterface tetgen (mesh); > > >> > > > >> > > > >> > // Determine the TRI3 elements which define the convex hull > > >> > tetgen.pointset_convexhull(); > > >> > > > >> > // find neighbors, etc > > >> > mesh.prepare_for_use(); > > >> > > > >> > > > >> > // Generate tetrahedra > > >> > tetgen.triangulate_conformingDelaunayMesh(X, Y); > > >> > > > >> > > > >> > > > >> > Now, here's where my trouble comes in: for any value of X and Y, > > >> > I get the following error: > > >> > > > >> > tetgen.C:9469: long int tetgenmesh::flip(tetgenmesh::queue*, > > >> > tetgenmesh::badface**): Assertion `fc != NONCONVEX' failed. > > >> > > > >> > > > >> > What am I doing wrong? > > >> > > > >> > Thanks, > > >> > John > > >> > > > >> > 
From: David Xu <dxu@my...>  20060628 22:21:03

Hi All, What are the definitions for X, Y values in tetgen.triangulate_conformingDelaunayMesh(X, Y); Thanks, David On 6/21/06, Steffen Petersen <steffen.petersen@...> wrote: > > > Find attached an example file. > > Steffen > > > David Xu schrieb: > > Hi Steffen and John, > > > > I've been trying to figure out to how to generate tetrahedral mesh in > > a simple cubic domain and I just happened to see your emails. I don't > > know much about tetgen. Could you guys send me your files so that I > > can learn from them? > > > > Thanks a lot! > > > > David > > > > On 6/8/06, Steffen Petersen <steffen.petersen@...> wrote: > > > >> > >> I cannot see anything suspicious in the code. > >> I have just tested a simliar example which worked for me. > >> I will send the file to you. > >> > >> However, I have noticed some other problems. > >> In optimized mode I had to change the CXXFLAGS > >> in contrib/tetgen/Makefile > >> > >> #CXXFLAGS += DTRILIBRARY > >> CXXFLAGS := DNDEBUG DTRILIBRARY fPIC > >> > >> in order to get the example working. > >> > >> Steffen > >> > >> > >> > >> > Hi, > >> > > >> > I'm trying to use tetgen to generate a tetrahedral > >> > mesh in the cube [1,1]^3. > >> > > >> > So far, my code is the following: > >> > > >> > > >> > Mesh mesh(3); > >> > > >> > // add nodes at the eight vertices of the cube > >> > ... > >> > > >> > // Instantiate TetGen interface > >> > TetGenMeshInterface tetgen (mesh); > >> > > >> > > >> > // Determine the TRI3 elements which define the convex hull > >> > tetgen.pointset_convexhull(); > >> > > >> > // find neighbors, etc > >> > mesh.prepare_for_use(); > >> > > >> > > >> > // Generate tetrahedra > >> > tetgen.triangulate_conformingDelaunayMesh(X, Y); > >> > > >> > > >> > > >> > Now, here's where my trouble comes in: for any value of X and Y, > >> > I get the following error: > >> > > >> > tetgen.C:9469: long int tetgenmesh::flip(tetgenmesh::queue*, > >> > tetgenmesh::badface**): Assertion `fc != NONCONVEX' failed. > >> > > >> > > >> > What am I doing wrong? > >> > > >> > Thanks, > >> > John > >> > > >> > > >> > _______________________________________________ > >> > Libmeshdevel mailing list > >> > Libmeshdevel@... > >> > https://lists.sourceforge.net/lists/listinfo/libmeshdevel > >> > > >> > >> > >> > >> _______________________________________________ > >> Libmeshdevel mailing list > >> Libmeshdevel@... > >> https://lists.sourceforge.net/lists/listinfo/libmeshdevel > >> > > > > 
From: David Knezevic <david.knezevic@ba...>  20060627 08:02:42

>> Every boundary node should be on a boundary side of some element  by >> what other definition is it a boundary node? > > A triangle can have a node which is on the boundary but no side > which is on the boundary. However, penalty BC should handle this > case with no problem. Yep, I was getting confused, I thought I had to do something special/ different with triangles because non boundary triangles can have a node on the boundary, but of course I don't because this node is also on a boundary side of another triangle. Thanks for pointing out my misconception. It turns out the weird behaviour I was seeing was caused because I was imposing of the _wrong_ boundary conditions for the PDE I was solving, not because of the way I was imposing the BCs. Thanks, David 
From: John Peterson <peterson@cf...>  20060627 01:25:17

Roy Stogner writes: > On Mon, 26 Jun 2006, David Knezevic wrote: > > > I was wondering what the preferred method is for imposing Dirichlet > > BCs on a domain with triangles (TRI6s in my case), > > My preferred method for setting Dirichlet BCs is by actually > integrating the penalized terms over each boundary side. Trying to > set BCs node by node on LAGRANGE elements is faster, but trying to do > it node by node on any other element is wrong. > > > since we have boundary nodes that are not on a "boundary side" of an > > element. > > Every boundary node should be on a boundary side of some element  by > what other definition is it a boundary node? A triangle can have a node which is on the boundary but no side which is on the boundary. However, penalty BC should handle this case with no problem. > Even if you're adding a penalty matrix node by node, you don't have to > add the penalty on every element containing each boundary node, you > just have to add it on some element containing each boundary node. If > a boundary node is part of three elements but it's only on a boundary > edge of two of them, that's fine  it doesn't matter if the final > system contains 2*penalty on the matrix diagonal and rhs or if it > contains 3*penalty; as long as the same value is on each, you'll still > get the right boundary value to within O(1/penalty) error. > > > Currently I'm looping over the elements and testing the location of > > each node to see if a BC should be imposed using the penalty method, > > i.e. I use the following code within the element loop: > > > > > > for(unsigned int n=0; n<elem>n_nodes(); n++) > > { > > const Point& p = elem>point(n); > > if(fabs(p(1)x2_max) < TOL) > > { > > Kuu(n,n) += penalty; // set u = 0 on this boundary > > Kvv(n,n) += penalty; > > Fv(n) += penalty; // set v = 1 on this boundary > > } > > } Do you definitely want to test p(1) against an "x2" value? I'm not sure at the numbering scheme is exactly. > > > > I think this should work (?), but it's giving me weird behaviour > > where the BC value seems to depend on the number of times the > > particular node is looped over. > > Could you be more specific? What BC values are you seeing along > y==x2_max? > > > I was wondering if there's a better way of doing this, and if there's > > anything wrong with the approach I described above. > > The approach above will break on nonLAGRANGE elements and on any > elements which use first order bases on a second order geometric > element. It assumes that every node at y==x2_max is a boundary node, > and it leaves the natural boundary conditions on every other boundary. > > None of those statements is necessarily a problem, though, so I'm not > sure what the bug could be. If by some chance, you do hit the same node with both a zero and a one boundary condition, then that node will get a value of 1/2 assigned to it. >  > Roy 
From: Roy Stogner <roystgnr@ic...>  20060626 21:20:49

On Mon, 26 Jun 2006, David Knezevic wrote: > I was wondering what the preferred method is for imposing Dirichlet > BCs on a domain with triangles (TRI6s in my case), My preferred method for setting Dirichlet BCs is by actually integrating the penalized terms over each boundary side. Trying to set BCs node by node on LAGRANGE elements is faster, but trying to do it node by node on any other element is wrong. > since we have boundary nodes that are not on a "boundary side" of an > element. Every boundary node should be on a boundary side of some element  by what other definition is it a boundary node? Even if you're adding a penalty matrix node by node, you don't have to add the penalty on every element containing each boundary node, you just have to add it on some element containing each boundary node. If a boundary node is part of three elements but it's only on a boundary edge of two of them, that's fine  it doesn't matter if the final system contains 2*penalty on the matrix diagonal and rhs or if it contains 3*penalty; as long as the same value is on each, you'll still get the right boundary value to within O(1/penalty) error. > Currently I'm looping over the elements and testing the location of > each node to see if a BC should be imposed using the penalty method, > i.e. I use the following code within the element loop: > > > for(unsigned int n=0; n<elem>n_nodes(); n++) > { > const Point& p = elem>point(n); > if(fabs(p(1)x2_max) < TOL) > { > Kuu(n,n) += penalty; // set u = 0 on this boundary > Kvv(n,n) += penalty; > Fv(n) += penalty; // set v = 1 on this boundary > } > } > > > I think this should work (?), but it's giving me weird behaviour > where the BC value seems to depend on the number of times the > particular node is looped over. Could you be more specific? What BC values are you seeing along y==x2_max? > I was wondering if there's a better way of doing this, and if there's > anything wrong with the approach I described above. The approach above will break on nonLAGRANGE elements and on any elements which use first order bases on a second order geometric element. It assumes that every node at y==x2_max is a boundary node, and it leaves the natural boundary conditions on every other boundary. None of those statements is necessarily a problem, though, so I'm not sure what the bug could be.  Roy 
From: David Knezevic <david.knezevic@ba...>  20060626 20:36:09

Hi all, I was wondering what the preferred method is for imposing Dirichlet BCs on a domain with triangles (TRI6s in my case), since we have boundary nodes that are not on a "boundary side" of an element. Currently I'm looping over the elements and testing the location of each node to see if a BC should be imposed using the penalty method, i.e. I use the following code within the element loop: for(unsigned int n=0; n<elem>n_nodes(); n++) { const Point& p = elem>point(n); if(fabs(p(1)x2_max) < TOL) { Kuu(n,n) += penalty; // set u = 0 on this boundary Kvv(n,n) += penalty; Fv(n) += penalty; // set v = 1 on this boundary } } I think this should work (?), but it's giving me weird behaviour where the BC value seems to depend on the number of times the particular node is looped over. I was wondering if there's a better way of doing this, and if there's anything wrong with the approach I described above. Thanks, David 
From: Lawrence Barlow <odwoqrlwpt@st...>  20060624 08:16:13

From: David Xu <dxu@my...>  20060623 22:53:56

Hi Ben, Thanks for your reply. I kinda understand the concept of the penalty scheme. The part I'm not sure is why we use L2 projection of the boundary data? Say, I just want to set the solution u=0 at the boudaries of a cubic domain and I have stiffness matrix Ke and mass matrix Me, the generalized eigenvalue problem equation is Ke*u=lamda*Me*u, the RHS is Me, not a vector. How would you impose penalty on Ke and Me to achieve u=0 BC? Some pseudo code would be even more helpful and highly appreciated. Regards, David On 6/23/06, Kirk, Benjamin (JSCEG) <benjamin.kirk1@...> wrote: > What you want is the linear system to reduce to the L2 projection of the > boundary data for the boundary degrees of freedom. The way we tend to > do that is to add the L2 projection, scaled by some *huge* penalty > value, to the matrix and RHS. In floating point arithmetic the idea is > that you add these huge values to the matrix and the other values in the > matrix are no longer significant and essentially are 0. > > One reason for doing this is that you only need to assess boundary > conditions in elements that have a full DIM1 face on the boundary and > can ignore elements with lowerdimensional boundary traces (edges in 3D, > nodes in 2D/3D). > > I am not familiar with your formulation, but I am almost certain you > want the penalty value multiplying the Me() term and not the Ke() term. > Again, the notion is to add the L2 projection of the boundary data, and > the mass matrix contains the terms necessary for the L2 projection on > the matrix side. > > Ben > > > > > Original Message > From: libmeshusersbounces@... > [mailto:libmeshusersbounces@...] On Behalf Of David > Xu > Sent: Thursday, June 22, 2006 4:59 PM > To: libmeshusers@... > Subject: [Libmeshusers] Question about Dirichlet BC imposed by > penaltymethod > > Hi, > > I know there are already many examples at libmesh web site using penalty > method to impose Dirichlet BC. However none of of them is generalized > eigenvalue problem related and most of them have vector as the RHS. My > problem is to solve the eigenvalue of a stationary 3D Schrodinger > equation. > > The element loop for the mass and stiffness matrices is the following: > >  >  > for (unsigned int qp=0; qp<qrule.n_points(); qp++) > { > const Real x = q_point[qp](0); > const Real y = q_point[qp](1); > const Real z = q_point[qp](2); > for (unsigned int i=0; i<phi.size(); i++) > for (unsigned int j=0; j<phi.size(); j++) > { > Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]; > Ke(i,j) += JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] + > (.5*pow(x, 2.)+.72*pow(y, 2.)+.845*pow(z, 2.))*phi[i][qp]*phi[j][qp]); > } > } >  >  > > How do I impose the simple Dirichlet BC( u=0 at the boundary) using > penalty method? Does the following code look right? > >  >  > for (unsigned int side=0; side<elem>n_sides(); side++) > if (elem>neighbor(side) == NULL) > { > const std::vector<std::vector<Real> >& phi_face = > fe_face>get_phi(); > const std::vector<std::vector<RealGradient> >& > dphi_face = fe_face>get_dphi(); > const std::vector<Real>& JxW_face = fe_face>get_JxW(); > const std::vector<Point >& qface_point = > fe_face>get_xyz(); > fe_face>reinit(elem, side); > > for (unsigned int qp=0; qp<qface.n_points(); qp++) > { > const Real xf = qface_point[qp](0); > const Real yf = qface_point[qp](1); > const Real zf = qface_point[qp](2); > > const Real penalty = 1.e10; > > for (unsigned int i=0; i<phi_face.size(); i++) > for (unsigned int j=0; > j<phi_face.size(); j++) > { > Me(i,j) += > JxW_face[qp]*0*penalty*phi_face[i][qp]*phi_face[j][qp]; > Ke(i,j) += > JxW_face[qp]*penalty*(.5*dphi_face[i][qp]*dphi_face[j][qp] + (.5*pow(xf, > 2.)+.72*pow(yf, 2.)+.845*pow(zf, 2.))*phi_face[i][qp]*phi_face[j][qp]); > } > } > } >  >  > > Thanks a lot in advance! > > David > > Using Tomcat but need to do more? Need to support web services, > security? > Get stuff done quickly with preintegrated technology to make your job > easier Download IBM WebSphere Application Server v.1.0.1 based on Apache > Geronimo > http://sel.asus.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642 > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 
From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20060623 17:18:40

What you want is the linear system to reduce to the L2 projection of the boundary data for the boundary degrees of freedom. The way we tend to do that is to add the L2 projection, scaled by some *huge* penalty value, to the matrix and RHS. In floating point arithmetic the idea is that you add these huge values to the matrix and the other values in the matrix are no longer significant and essentially are 0. One reason for doing this is that you only need to assess boundary conditions in elements that have a full DIM1 face on the boundary and can ignore elements with lowerdimensional boundary traces (edges in 3D, nodes in 2D/3D). =20 I am not familiar with your formulation, but I am almost certain you want the penalty value multiplying the Me() term and not the Ke() term. Again, the notion is to add the L2 projection of the boundary data, and the mass matrix contains the terms necessary for the L2 projection on the matrix side. Ben =20 Original Message From: libmeshusersbounces@... [mailto:libmeshusersbounces@...] On Behalf Of David Xu Sent: Thursday, June 22, 2006 4:59 PM To: libmeshusers@... Subject: [Libmeshusers] Question about Dirichlet BC imposed by penaltymethod Hi, I know there are already many examples at libmesh web site using penalty method to impose Dirichlet BC. However none of of them is generalized eigenvalue problem related and most of them have vector as the RHS. My problem is to solve the eigenvalue of a stationary 3D Schrodinger equation. The element loop for the mass and stiffness matrices is the following:   for (unsigned int qp=3D0; qp<qrule.n_points(); qp++) { const Real x =3D q_point[qp](0); const Real y =3D q_point[qp](1); const Real z =3D q_point[qp](2); for (unsigned int i=3D0; i<phi.size(); i++) for (unsigned int j=3D0; j<phi.size(); j++) { Me(i,j) +=3D JxW[qp]*phi[i][qp]*phi[j][qp]; Ke(i,j) +=3D JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] + (.5*pow(x, 2.)+.72*pow(y, 2.)+.845*pow(z, 2.))*phi[i][qp]*phi[j][qp]); } }   How do I impose the simple Dirichlet BC( u=3D0 at the boundary) using penalty method? Does the following code look right?   for (unsigned int side=3D0; side<elem>n_sides(); side++) if (elem>neighbor(side) =3D=3D NULL) { const std::vector<std::vector<Real> >& phi_face =3D fe_face>get_phi(); const std::vector<std::vector<RealGradient> >& dphi_face =3D fe_face>get_dphi(); const std::vector<Real>& JxW_face =3D = fe_face>get_JxW(); const std::vector<Point >& qface_point =3D fe_face>get_xyz(); fe_face>reinit(elem, side); for (unsigned int qp=3D0; qp<qface.n_points(); qp++) { const Real xf =3D qface_point[qp](0); const Real yf =3D qface_point[qp](1); const Real zf =3D qface_point[qp](2); const Real penalty =3D 1.e10; for (unsigned int i=3D0; i<phi_face.size(); = i++) for (unsigned int j=3D0; j<phi_face.size(); j++) { Me(i,j) +=3D JxW_face[qp]*0*penalty*phi_face[i][qp]*phi_face[j][qp]; Ke(i,j) +=3D JxW_face[qp]*penalty*(.5*dphi_face[i][qp]*dphi_face[j][qp] + (.5*pow(xf, 2.)+.72*pow(yf, 2.)+.845*pow(zf, 2.))*phi_face[i][qp]*phi_face[j][qp]); } } }   Thanks a lot in advance! David Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with preintegrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.asus.falkag.net/sel?cmd=3Dlnk&kid=3D120709&bid=3D263057&dat=3D= 121642 _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: David Xu <dxu@my...>  20060622 21:58:52

Hi, I know there are already many examples at libmesh web site using penalty method to impose Dirichlet BC. However none of of them is generalized eigenvalue problem related and most of them have vector as the RHS. My problem is to solve the eigenvalue of a stationary 3D Schrodinger equation. The element loop for the mass and stiffness matrices is the following:  for (unsigned int qp=0; qp<qrule.n_points(); qp++) { const Real x = q_point[qp](0); const Real y = q_point[qp](1); const Real z = q_point[qp](2); for (unsigned int i=0; i<phi.size(); i++) for (unsigned int j=0; j<phi.size(); j++) { Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]; Ke(i,j) += JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] + (.5*pow(x, 2.)+.72*pow(y, 2.)+.845*pow(z, 2.))*phi[i][qp]*phi[j][qp]); } }  How do I impose the simple Dirichlet BC( u=0 at the boundary) using penalty method? Does the following code look right?  for (unsigned int side=0; side<elem>n_sides(); side++) if (elem>neighbor(side) == NULL) { const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi(); const std::vector<std::vector<RealGradient> >& dphi_face = fe_face>get_dphi(); const std::vector<Real>& JxW_face = fe_face>get_JxW(); const std::vector<Point >& qface_point = fe_face>get_xyz(); fe_face>reinit(elem, side); for (unsigned int qp=0; qp<qface.n_points(); qp++) { const Real xf = qface_point[qp](0); const Real yf = qface_point[qp](1); const Real zf = qface_point[qp](2); const Real penalty = 1.e10; for (unsigned int i=0; i<phi_face.size(); i++) for (unsigned int j=0; j<phi_face.size(); j++) { Me(i,j) += JxW_face[qp]*0*penalty*phi_face[i][qp]*phi_face[j][qp]; Ke(i,j) += JxW_face[qp]*penalty*(.5*dphi_face[i][qp]*dphi_face[j][qp] + (.5*pow(xf, 2.)+.72*pow(yf, 2.)+.845*pow(zf, 2.))*phi_face[i][qp]*phi_face[j][qp]); } } }  Thanks a lot in advance! David 
From: Roy Stogner <roystgnr@ic...>  20060621 16:36:17

On Wed, 21 Jun 2006, li pan wrote: > just one simple question. Does it make any difference, > if solution_in was or was not set to zero before I > call PetscLinearSolver::solve(...)? It should. The solution you pass in gets used as the initial iterate to the solver, and so should be as close to the expected final solution as possible. That may mean you want to set it to zero (e.g. if you're solving for the difference between two Newton iterates), but it may mean you want to leave it set to a previous iterate (e.g. if you're doing one linear solve for the solution at the next timestep, or if you're using a nonlinear iteration scheme and solving directly for the next iterate).  Roy 
From: li pan <li76pan@ya...>  20060621 16:28:04

Dear developers, just one simple question. Does it make any difference, if solution_in was or was not set to zero before I call PetscLinearSolver::solve(...)? regards pan __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com 
From: li pan <li76pan@ya...>  20060620 08:41:20

thanx Roy, it works with static test. But I still have problem with dynamic test, because there are more matrices and vectors (stiffness, mass, force, displacement, velocity and acceleration). The true system.matrix and system.rhs are all effective one. They will be calculated after assembly. Only mass, stiffness and force are assembled. So I called dof_map.constrain_XXXX for those three and got the following error. Hereby, one thing seems odd for is both stiffness and mass have value at diagonal number 3035, but system.matrix doesn't. system.matrix=stiffness+A0*mass (A0 is a constant scale coefficient) [0]PETSC ERROR: MatMissingDiagonal_SeqAIJ() line 998 in src/mat/impls/aij/seq/aij.c [0]PETSC ERROR: Petsc has generated inconsistent data! [0]PETSC ERROR: Matrix is missing diagonal number 3035! [0]PETSC ERROR: MatILUFactorSymbolic_SeqAIJ() line 886 in src/mat/impls/aij/seq/aijfact.c [0]PETSC ERROR: MatILUFactorSymbolic() line 4312 in src/mat/interface/matrix.c [0]PETSC ERROR: PCSetUp_ILU() line 779 in src/ksp/pc/impls/factor/ilu/ilu.c [0]PETSC ERROR: PCSetUp() line 798 in src/ksp/pc/interface/precon.c [0]PETSC ERROR: KSPSetUp() line 234 in src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: KSPSolve() line 332 in src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: User provided function() line 288 in unknowndirectory/src/numerics/petsc_linear_solver.C On Mon, 19 Jun 2006, li pan wrote: > I used error estimator and mesh refinement for the > first time. I wonder how to define the three > attributes, namely, refine_fraction, coarsen_fraction, > and max_level. I think for special problem, the value > must be different. I did some test. Tried with > different combination. The Newton iteration which > converged without mesh refinement well before, can't > converge with some refinement. What could likely be > the reason? Well, the reason that just killed one of our visitors to UT Austin is simple: In the matrix assembly stage, after you've built the local dense matrix and vector for a given cell, but before you add them to the global sparse matrix and vector, you need to call something like: dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices); We really need to get some code into libMesh to warn users about this potential error, but I'm not sure how to do it. Maybe a static boolean in DofMap that initializes to false, gets turned on by the constrain_* member functions, and then gets tested by the System::solve functions whenever DofMap::n_constrained_dofs() > 0? The other theoretical problem with a Newton iteration is that if you're solving for the delta between Newton iterates, it may be possible that an inexact linear solve will give you illconstrained hanging DoFs and that subsequent Newton steps will be unable to correct them. I call this a "theoretical problem" because I've never seen it in practice, but I think I may code a DofMap::enforce_constraints() method (to be called by System::solve after the linear solver is done) to be safe in the future.  Roy Stogner __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com 
From: Roy Stogner <roystgnr@ic...>  20060619 16:10:49

On Mon, 19 Jun 2006, li pan wrote: > I used error estimator and mesh refinement for the > first time. I wonder how to define the three > attributes, namely, refine_fraction, coarsen_fraction, > and max_level. I think for special problem, the value > must be different. I did some test. Tried with > different combination. The Newton iteration which > converged without mesh refinement well before, can't > converge with some refinement. What could likely be > the reason? Well, the reason that just killed one of our visitors to UT Austin is simple: In the matrix assembly stage, after you've built the local dense matrix and vector for a given cell, but before you add them to the global sparse matrix and vector, you need to call something like: dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices); We really need to get some code into libMesh to warn users about this potential error, but I'm not sure how to do it. Maybe a static boolean in DofMap that initializes to false, gets turned on by the constrain_* member functions, and then gets tested by the System::solve functions whenever DofMap::n_constrained_dofs() > 0? The other theoretical problem with a Newton iteration is that if you're solving for the delta between Newton iterates, it may be possible that an inexact linear solve will give you illconstrained hanging DoFs and that subsequent Newton steps will be unable to correct them. I call this a "theoretical problem" because I've never seen it in practice, but I think I may code a DofMap::enforce_constraints() method (to be called by System::solve after the linear solver is done) to be safe in the future.  Roy Stogner 
From: Roy Stogner <roystgnr@ic...>  20060619 16:02:29

On Mon, 19 Jun 2006, Guillaume ANCIAUX wrote: > I was wondering is there is a way to obtain a global numbering of the nodes > that does not change with partitions. The node numbering should stay the same unless you remove nodes, shouldn't it? I know we've got more than a little code in some of the finite element classes that assumes node ordering is never changed. Keep in mind that the node numbering and the degree of freedom numbering aren't the same thing. I don't think there's a natural way to keep the degree of freedom numbering constant through a repartitioning, since we want to keep a contiguous range of degree of freedom indices on each processor, and trying to repartition while maintaining that state and without renumbering would be too limiting. Can you create such a numbering yourself? If you're using isoparametric Lagrange variables, then add another one, and assign to each of its indices a unique value. After any repartitioning, the System::project_vector calls should make sure that your lookup variable gets updated correctly.  Roy Stogner 