From: Joa L. <li...@jo...> - 2009-10-13 13:17:26
|
Hi all, I would like to use libmesh to solve a rather simple equation (Laplace/poisson) but in a domain with a somewhat funny shape of the boundary. To do this I created a mesh using gmsh, modified example 14 a bit so it reads my mesh instead of the l-shaped domain. My problem is that when I refine I get "flat" surfaces that does not follow my boundary (which is not flat;). How do I go about and move the new boundary vertices of my tets (I use tets) to the real boundary? As for the geometry etc. I now how to do it, I`m just not so familiar with libmesh. kind regards Joa Ljungvall |
From: David K. <dknez@MIT.EDU> - 2009-10-13 13:25:40
|
I think the simplest thing to do is use second order elements in gmsh. Then when you refine the new nodes will interpolate the quadratic approximation to the boundary of your domain. - Dave Joa Ljungvall wrote: > Hi all, > > I would like to use libmesh to solve a rather simple equation (Laplace/poisson) > but in a domain with a somewhat funny shape of the boundary. To do this I > created a mesh using gmsh, modified example 14 a bit so it reads my mesh > instead of the l-shaped domain. My problem is that when I refine I get "flat" > surfaces that does not follow my boundary (which is not flat;). How do I go > about and move the new boundary vertices of my tets (I use tets) to the real > boundary? As for the geometry etc. I now how to do it, I`m just not so familiar > with libmesh. > > > kind regards > > Joa Ljungvall > > > ------------------------------------------------------------------------------ > Come build with us! The BlackBerry(R) Developer Conference in SF, CA > is the only developer event you need to attend this year. Jumpstart your > developing skills, take BlackBerry mobile applications to market and stay > ahead of the curve. Join us from November 9 - 12, 2009. Register now! > http://p.sf.net/sfu/devconference > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > |
From: Joa L. <li...@jo...> - 2009-10-13 14:02:00
|
Hi, I'm not sure I understood the answer. I want to start with a very course mesh, and refine it based on the solution on the course(r) mesh, i.e. adaptive refinment. So using gmsh I would have to communicate the solution in each step and tell gmsh which regions to refine, or? And if so, how? Further more, my domain have some regions where a cylinder it cut by a hexagonal, so a quadratic approximation would not do so well, would it? cheers Joa On Tue, Oct 13, 2009 at 09:24:53AM -0400, David Knezevic wrote: > I think the simplest thing to do is use second order elements in gmsh. Then > when you refine the new nodes will interpolate the quadratic > approximation to > the boundary of your domain. > > - Dave > > > > Joa Ljungvall wrote: > >Hi all, > > > >I would like to use libmesh to solve a rather simple equation (Laplace/poisson) > >but in a domain with a somewhat funny shape of the boundary. To do > >this I created a mesh using gmsh, modified example 14 a bit so it > >reads my mesh > >instead of the l-shaped domain. My problem is that when I refine I get "flat" > >surfaces that does not follow my boundary (which is not flat;). How do I go > >about and move the new boundary vertices of my tets (I use tets) to the real > >boundary? As for the geometry etc. I now how to do it, I`m just not so familiar > >with libmesh. > > > > > >kind regards > > > >Joa Ljungvall > > > > > >------------------------------------------------------------------------------ > >Come build with us! The BlackBerry(R) Developer Conference in SF, CA > >is the only developer event you need to attend this year. Jumpstart your > >developing skills, take BlackBerry mobile applications to market > >and stay ahead of the curve. Join us from November 9 - 12, 2009. > >Register now! > >http://p.sf.net/sfu/devconference > >_______________________________________________ > >Libmesh-users mailing list > >Lib...@li... > >https://lists.sourceforge.net/lists/listinfo/libmesh-users |
From: David K. <dknez@MIT.EDU> - 2009-10-13 13:59:35
|
If you generate a mesh in gmsh that is sufficiently fine that a mesh of second order elements captures the geometry well enough, then you can read that mesh into libMesh, and just use libMesh's adaptive refinement as in the examples. In that situation, libMesh's adaptive refinement will interpolate the quadratic boundaries when you refine the boundary elements... But if you want to start with a very coarse mesh, then you might have to move boundary nodes around etc. I've never had to do that myself, so perhaps someone else on the list can offer you some choice advice. - Dave Joa Ljungvall wrote: > Hi, > > I'm not sure I understood the answer. I want to start with a very course mesh, > and refine it based on the solution on the course(r) mesh, i.e. adaptive > refinment. So using gmsh I would have to communicate the solution in each step > and tell gmsh which regions to refine, or? And if so, how? Further more, my > domain have some regions where a cylinder it cut by a hexagonal, so a > quadratic approximation would not do so well, would it? > > > cheers > > > Joa > > On Tue, Oct 13, 2009 at 09:24:53AM -0400, David Knezevic wrote: > >> I think the simplest thing to do is use second order elements in gmsh. Then >> when you refine the new nodes will interpolate the quadratic >> approximation to >> the boundary of your domain. >> >> - Dave >> >> >> >> Joa Ljungvall wrote: >> >>> Hi all, >>> >>> I would like to use libmesh to solve a rather simple equation (Laplace/poisson) >>> but in a domain with a somewhat funny shape of the boundary. To do >>> this I created a mesh using gmsh, modified example 14 a bit so it >>> reads my mesh >>> instead of the l-shaped domain. My problem is that when I refine I get "flat" >>> surfaces that does not follow my boundary (which is not flat;). How do I go >>> about and move the new boundary vertices of my tets (I use tets) to the real >>> boundary? As for the geometry etc. I now how to do it, I`m just not so familiar >>> with libmesh. >>> >>> >>> kind regards >>> >>> Joa Ljungvall >>> >>> >>> ------------------------------------------------------------------------------ >>> Come build with us! The BlackBerry(R) Developer Conference in SF, CA >>> is the only developer event you need to attend this year. Jumpstart your >>> developing skills, take BlackBerry mobile applications to market >>> and stay ahead of the curve. Join us from November 9 - 12, 2009. >>> Register now! >>> http://p.sf.net/sfu/devconference >>> _______________________________________________ >>> Libmesh-users mailing list >>> Lib...@li... >>> https://lists.sourceforge.net/lists/listinfo/libmesh-users >>> |
From: Joa L. <li...@jo...> - 2009-10-13 15:36:28
|
Hi, The problem is that I have to solve my equation on the same domain with very different boundary conditions, making a uniform refinement very inefficient and a waste of time, and worse in my case, RAM. So I do need to start with a very course mesh and refine. cheers Joa On Tue, Oct 13, 2009 at 09:55:54AM -0400, David Knezevic wrote: > If you generate a mesh in gmsh that is sufficiently fine that a mesh > of second order elements captures the geometry well enough, then you > can read that mesh into libMesh, and just use libMesh's adaptive > refinement as in the examples. In that situation, libMesh's adaptive > refinement will interpolate the quadratic boundaries when you refine > the boundary elements... > > But if you want to start with a very coarse mesh, then you might > have to move boundary nodes around etc. I've never had to do that > myself, so perhaps someone else on the list can offer you some > choice advice. > > - Dave > > > Joa Ljungvall wrote: > >Hi, > > > >I'm not sure I understood the answer. I want to start with a very > >course mesh, and refine it based on the solution on the course(r) > >mesh, i.e. adaptive > >refinment. So using gmsh I would have to communicate the solution > >in each step and tell gmsh which regions to refine, or? And if so, > >how? Further more, my domain have some regions where a cylinder it > >cut by a hexagonal, so a quadratic approximation would not do so > >well, would it? > > > > > >cheers > > > > > >Joa > > > >On Tue, Oct 13, 2009 at 09:24:53AM -0400, David Knezevic wrote: > >>I think the simplest thing to do is use second order elements in gmsh. Then > >>when you refine the new nodes will interpolate the quadratic > >>approximation to > >>the boundary of your domain. > >> > >>- Dave > >> > >> > >> > >>Joa Ljungvall wrote: > >>>Hi all, > >>> > >>>I would like to use libmesh to solve a rather simple equation (Laplace/poisson) > >>>but in a domain with a somewhat funny shape of the boundary. To do > >>>this I created a mesh using gmsh, modified example 14 a bit so it > >>>reads my mesh > >>>instead of the l-shaped domain. My problem is that when I refine I get "flat" > >>>surfaces that does not follow my boundary (which is not flat;). How do I go > >>>about and move the new boundary vertices of my tets (I use tets) to the real > >>>boundary? As for the geometry etc. I now how to do it, I`m just not so familiar > >>>with libmesh. > >>> > >>> > >>>kind regards > >>> > >>>Joa Ljungvall > >>> > >>> > >>>------------------------------------------------------------------------------ > >>>Come build with us! The BlackBerry(R) Developer Conference in SF, CA > >>>is the only developer event you need to attend this year. Jumpstart your > >>>developing skills, take BlackBerry mobile applications to market > >>>and stay ahead of the curve. Join us from November 9 - 12, 2009. > >>>Register now! > >>>http://p.sf.net/sfu/devconference > >>>_______________________________________________ > >>>Libmesh-users mailing list > >>>Lib...@li... > >>>https://lists.sourceforge.net/lists/listinfo/libmesh-users |
From: John P. <pet...@cf...> - 2009-10-13 15:55:32
|
On Tue, Oct 13, 2009 at 10:36 AM, Joa Ljungvall <li...@jo...> wrote: > Hi, > > The problem is that I have to solve my equation on the same domain with > very different boundary conditions, making a uniform refinement very > inefficient and a waste of time, and worse in my case, RAM. So I do need > to start with a very course mesh and refine. As David explained, there is unfortunately currently no "feedback" mechanism between the mesh refinement procedure and an underlying geometric description. One approach might be to maintain a suitably fine, lower-dimensional (and hence, low memory) discretization of the boundary as a secondary Mesh object in LibMesh, and then somehow use this fine discretization when placing new points, but again, there is no existing mechanism for doing this currently in the library. -- John |
From: Joa L. <li...@jo...> - 2009-10-13 16:41:07
|
Hi again, during the day I been working a bit on the problem: Currently I have the exact geometry inside my program (implemented via CGS and GEANT4). What I have tried is to, after the refinement move the new boundary points using a code like this one (it tries to move all nodes on a boundary to the correct "radius" along the rhat vector void libMESH::AdjustMesh(Mesh &mesh,BaseBoundary *BC) { MeshBase::element_iterator el = mesh.active_elements_begin(); const MeshBase::element_iterator end_el = mesh.active_elements_end(); for ( ; el != end_el; ++el) { for (unsigned int s=0; s<(*el)->n_sides(); s++) if ((*el)->neighbor(s) == NULL) { AutoPtr<DofObject> face = (*el)->side(s); for(int atnode=0; atnode<3; ++atnode){ Node *node = (dynamic_cast<Elem*>(face.get())->set_node(atnode)); //if z==0 its ok if(fabs((*node)(2))>1e-6){ if((*node)(2)<BC->CylLength()){//Simple case //Check which boundary is the closest double Router = BC->RAtZ(MakeHepVector((*node))); double Rinner = BC->RadiusOnCentralContact(MakeHepVector((*node))); double Rnode = sqrt((*node)(0)*(*node)(0)+(*node)(1)*(*node)(1)); if(fabs(Router-Rnode)<fabs(Rinner-Rnode)){ (*node)(0)*=Router/Rnode; (*node)(1)*=Router/Rnode; } else { (*node)(0)*=Rinner/Rnode; (*node)(1)*=Rinner/Rnode; } } else {//Difficult, we have to check if the node should be on //one of the "circles" } } } } } } where BaseBoundary is a class to take care of the geometry. Trying with this I get EquationSystems n_systems()=1 System "Laplace" Type "LinearImplicit" Variables="u" Finite Element Types="LAGRANGE" Approximation Orders="FIRST" n_dofs()=715 n_local_dofs()=715 n_constrained_dofs()=0 n_vectors()=1 Beginning Solve 0 System has: 715 degrees of freedom. Linear solver converged at step: 5000, final residual: 1.1091e-16 Refining the mesh... Beginning Solve 1 ERROR: negative Jacobian: -83.5253 in element 7069 [0] src/fe/fe_map.C, line 341, compiled Oct 11 2009 at 15:03:39 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Abort trap It is quite possible that I move the wrong nodes the wrong way, but is my problem ;). The general question is, is this a possible route to take? I call this routine after refinement and coursening, before equation_systems.reinit(). cheers Joa On Tue, Oct 13, 2009 at 10:54:57AM -0500, John Peterson wrote: > On Tue, Oct 13, 2009 at 10:36 AM, Joa Ljungvall <li...@jo...> wrote: > > Hi, > > > > The problem is that I have to solve my equation on the same domain with > > very different boundary conditions, making a uniform refinement very > > inefficient and a waste of time, and worse in my case, RAM. So I do need > > to start with a very course mesh and refine. > > As David explained, there is unfortunately currently no "feedback" > mechanism between the mesh refinement procedure and an underlying > geometric description. > > One approach might be to maintain a suitably fine, lower-dimensional > (and hence, low memory) discretization of the boundary as a secondary > Mesh object in LibMesh, and then somehow use this fine discretization > when placing new points, but again, there is no existing mechanism for > doing this currently in the library. > > -- > John |
From: John P. <pet...@cf...> - 2009-10-13 16:59:06
|
On Tue, Oct 13, 2009 at 11:40 AM, Joa Ljungvall <li...@jo...> wrote: > It is quite possible that I move the wrong nodes the wrong way, but is my > problem ;). The general question is, is this a possible route to take? I call > this routine after refinement and coursening, before equation_systems.reinit(). I seem to remember there being some restriction when doing refinement/coarsening and mesh redistribution, i.e. all mesh movement had to be performed before any refinement/coarsening. I think it had to do with the coarsening step, but the way we handle constraints also assumes some level of parent is conforming with the neighbor... and I don't think this is the case with mesh movement. -- John |
From: Roy S. <roy...@ic...> - 2009-10-13 17:06:06
|
On Tue, 13 Oct 2009, John Peterson wrote: > On Tue, Oct 13, 2009 at 11:40 AM, Joa Ljungvall <li...@jo...> wrote: >> It is quite possible that I move the wrong nodes the wrong way, but is my >> problem ;). The general question is, is this a possible route to take? I call >> this routine after refinement and coursening, before equation_systems.reinit(). > > I seem to remember there being some restriction when doing > refinement/coarsening and mesh redistribution, i.e. all mesh movement > had to be performed before any refinement/coarsening. I think it had > to do with the coarsening step, but the way we handle constraints also > assumes some level of parent is conforming with the neighbor... and I > don't think this is the case with mesh movement. The parent can always be made conforming with the neighbor in the case of a hanging node constraint - you just need to note that the hanging node's location is constrained as well, and keep that constraint satisfied. But failure to do that would just lead to potentially inaccurate projections and constraint equations. I'm not sure how someone would run into an inverted Jacobian... unless the coarse mesh was so far from the real boundary that snapping nodes to the boundary resulted in a completely twisted-around element. --- Roy |
From: Joa L. <li...@jo...> - 2009-10-13 17:14:18
|
Hi again, If I: 0) Solve the equations one time... 1) Refine my mesh 2) Move my nodes 3) Copies all active elements to a new mesh, i.e. my new mesh have the right nodes locations etc. but no "refinement" (is this true if I do it right?) 4) I change mesh in my equation system 5) build my matrix and RHS, solves 6) -> 1 until pleased It is clearly time consuming (alloc free system calls) and it will require almost twice the memory... apart from that, could it work? cheers Joa On Tue, Oct 13, 2009 at 11:58:30AM -0500, John Peterson wrote: > On Tue, Oct 13, 2009 at 11:40 AM, Joa Ljungvall <li...@jo...> wrote: > > It is quite possible that I move the wrong nodes the wrong way, but is my > > problem ;). The general question is, is this a possible route to take? I call > > this routine after refinement and coursening, before equation_systems.reinit(). > > I seem to remember there being some restriction when doing > refinement/coarsening and mesh redistribution, i.e. all mesh movement > had to be performed before any refinement/coarsening. I think it had > to do with the coarsening step, but the way we handle constraints also > assumes some level of parent is conforming with the neighbor... and I > don't think this is the case with mesh movement. > > -- > John |
From: Roy S. <roy...@ic...> - 2009-10-13 17:17:13
|
> 3) Copies all active elements to a new mesh, i.e. my new mesh have the right > nodes locations etc. but no "refinement" (is this true if I do it right?) This won't work if there's any adaptive refinement in the original mesh - libMesh relies on the parent element hierarchy to figure out hanging node constraints. --- Roy |
From: Joa L. <li...@jo...> - 2009-10-13 17:36:21
|
Hmmm... If I change the locations of the new nodes when they are made, i.e. somewhere deep down in the MeshRefinement class, can I do what I would like to do? If I've understood everything correct a tet is subdivided by adding vertices halfway between the old ones. But there is nothing magic about that, or? One could choose to put the vertices somewhere else, as long as the element remains valid. Coarsening is of course much simpler and can remain the way it is. If this could work in princple, I'm willing to give it a try. A problem would be to define an interface with a geometry, but if I get to that... cheers Joa On Tue, Oct 13, 2009 at 12:17:03PM -0500, Roy Stogner wrote: > > >3) Copies all active elements to a new mesh, i.e. my new mesh have the right > > nodes locations etc. but no "refinement" (is this true if I do it right?) > > This won't work if there's any adaptive refinement in the original > mesh - libMesh relies on the parent element hierarchy to figure out > hanging node constraints. > --- > Roy |
From: Roy S. <roy...@ic...> - 2009-10-13 17:47:41
|
On Tue, 13 Oct 2009, Joa Ljungvall wrote: > If I change the locations of the new nodes when they are made, i.e. somewhere > deep down in the MeshRefinement class, can I do what I would like to do? I don't think that would make any difference. If moving nodes eventually gives you an inverted element, it probably won't matter exactly when the nodes are moved. I do think it would be a big mistake to start planning major low-level changes before we fully understand why the less intrusive algorithm is failing. Maybe get at that negative Jacobian failure in the debugger: is the inverted element active? A parent? If the latter, what code is calculating the Jacobian, and can it be worked around somehow? --- Roy |
From: Joa L. <li...@jo...> - 2009-10-13 18:15:22
|
Hi again, Ok. Seems sensible to me. Now, I'm a nuclear physicist so FEM is not quite what I do for a living, but presently I'm giving it a try to see if it can improve some detector simulations we are doing, presently done with finite difference solvers, so bare with me even if I'm not 100% on target all the time. I`ve already tried with dealII and that worked ok, but for good or bad reasons (I'm not sure;) I want simple tets with first order shape functions, hence libmesh. As for getting inverted element, this means I've moved a node to the other side of the opposite face of my tet? This could very well be a bug in my code. I can, and will, check that I don't move nodes in an unreasnoble way. Also, I'll get back with a stack/call trace and the nodes of the failing element etc. But Paris time is 20:13 and I got other duties calling for me so I guess that will be tomorrows project. cheers Joa On Tue, Oct 13, 2009 at 12:47:26PM -0500, Roy Stogner wrote: > > On Tue, 13 Oct 2009, Joa Ljungvall wrote: > > >If I change the locations of the new nodes when they are made, i.e. somewhere > >deep down in the MeshRefinement class, can I do what I would like to do? > > I don't think that would make any difference. If moving nodes > eventually gives you an inverted element, it probably won't matter > exactly when the nodes are moved. > > I do think it would be a big mistake to start planning major low-level > changes before we fully understand why the less intrusive algorithm is > failing. Maybe get at that negative Jacobian failure in the debugger: > is the inverted element active? A parent? If the latter, what code > is calculating the Jacobian, and can it be worked around somehow? > --- > Roy |
From: Roy S. <roy...@ic...> - 2009-10-13 18:19:11
|
On Tue, 13 Oct 2009, Joa Ljungvall wrote: > As for getting inverted element, this means I've moved a node to the other > side of the opposite face of my tet? This could very well be a bug in my > code. I can, and will, check that I don't move nodes in an unreasnoble way. If you're sticking to first order shape functions, make sure you use first order geometric elements as well. Moving a node past the opposite face is the only way to invert a first-order tet, but things get more complicated when the mapping function is quadratic. --- Roy |
From: Joa L. <li...@jo...> - 2009-10-13 21:49:47
|
The problem with a negative jacobian seems solved. It was indeed a node that was moved in the wrong way by my code. Now the program iterates two times but then I get (running in gdb) EquationSystems n_systems()=1 System "Laplace" Type "LinearImplicit" Variables="u" Finite Element Types="LAGRANGE" Approximation Orders="FIRST" n_dofs()=715 n_local_dofs()=715 n_constrained_dofs()=0 n_vectors()=1 Beginning Solve 0 System has: 715 degrees of freedom. Linear solver converged at step: 5000, final residual: 1.1091e-16 Refining the mesh... Beginning Solve 1 System has: 1043 degrees of freedom. Linear solver converged at step: 5000, final residual: 9.00376e-17 Refining the mesh... [0] /usr/local/libmesh/include/geom/remote_elem.h, line 91, compiled Oct 11 2009 at 15:02:42 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Program received signal SIGABRT, Aborted. 0x00007fff80dcdff6 in __kill () (gdb) bt #0 0x00007fff80dcdff6 in __kill () #1 0x00007fff80e6f072 in abort () #2 0x00007fff811185d2 in __gnu_cxx::__verbose_terminate_handler () #3 0x00007fff81116ae1 in __cxxabiv1::__terminate () #4 0x00007fff81116b16 in std::terminate () #5 0x00007fff81116bfc in __cxa_throw () #6 0x000000010512afa4 in RemoteElem::n_nodes () #7 0x000000010510acf4 in DofMap::dof_indices () #8 0x000000010510e2d9 in DofMap::add_neighbors_to_send_list () #9 0x000000010511122f in DofMap::distribute_dofs () #10 0x000000010548509d in EquationSystems::reinit () #11 0x0000000101d65d5a in AGATAGeFEM::libMESH::LaplaceProblem::run (this=0x1019061a0) at ../../libs/libmesh/LaplaceProblem.cc:554 #12 0x0000000100004e67 in DoFields (deffile=<value temporarily unavailable, due to optimizations>, solidfile=<value temporarily unavailable, due to optimizations>, OutputBase=0x7fff5fbff6ad "test", gnuplot=false, refine=0.90000000000000002, course=0.10000000000000001) at main.cc:454 #13 0x0000000100008afd in operator==<char, std::char_traits<char>, std::allocator<char> > [inlined] () at /usr/include/c++/4.2.1/bits/basic_string.h:1247 #14 0x0000000100008afd in main (argc=5, argv=0x7fff5fbff510) at main.cc:1249 This I don't understand at all... Joa On Tue, Oct 13, 2009 at 01:18:49PM -0500, Roy Stogner wrote: > > > On Tue, 13 Oct 2009, Joa Ljungvall wrote: > > >As for getting inverted element, this means I've moved a node to the other > >side of the opposite face of my tet? This could very well be a bug in my > >code. I can, and will, check that I don't move nodes in an unreasnoble way. > > If you're sticking to first order shape functions, make sure you use > first order geometric elements as well. Moving a node past the > opposite face is the only way to invert a first-order tet, but things > get more complicated when the mapping function is quadratic. > --- > Roy |
From: Roy S. <roy...@ic...> - 2009-10-13 22:09:19
|
On Tue, 13 Oct 2009, Joa Ljungvall wrote: > #5 0x00007fff81116bfc in __cxa_throw () > #6 0x000000010512afa4 in RemoteElem::n_nodes () > #7 0x000000010510acf4 in DofMap::dof_indices () > #8 0x000000010510e2d9 in DofMap::add_neighbors_to_send_list () > > This I don't understand at all... Hmm... there was a bug last year that looked like this, but it was fixed before 0.6.3 was released. It may be that you've inadvertently triggered the same problem that the bug did, though: When we refine an element, we need to create any new nodes its children will have, or we need to point them to the appropriate existing node in the mesh. We look for that existing node by location: calculate the point where we'd build a new node, see if any node exists at that point yet, if so attach to it and if not then build a new one. But what happens if you refine one element, create a hanging node, then *move the hanging node*? When you refine the element's neighbor, it will look to see if a node exists yet on that edge or face, the answer will be "no", and a new node will be created. This tears the topology of the mesh, which is a nasty bug that can manifest as the incorrect creation of a RemoteElem in our find_neighbors() algorithm. (but which should have triggered a libmesh_assert() earlier - did you run in devel or dbg mode?) Anyway, the solution is probably to, after you've looped over all nodes, loop over all constraint equations in the dof map and apply them to nodes' xyz coordinates. --- Roy |
From: Roy S. <roy...@ic...> - 2009-10-26 17:21:50
|
On Mon, 26 Oct 2009, Joa Ljungvall wrote: > What I'm trying to do is: > > 1) Refine the mesh > 2) Find and store hanging nodes > 3) Move points on boundary to the "real" boundary > then I do a do{...}While() including > 4) Check that Jacobian>0, if not switch node 0 and 2 (the pointers in the elem) > I'm not sure this doesn't mess up something for neighbors.. This merely hides a symptom without fixing the problem. If you have two elements ABC and CBD, and you invert one of them by moving a node too far: A---------C A---------C \ /| \ .-'/ \ / | \ D / \ / | \ | / \ / | \|/ B----D B Changing CBD into CDB doesn't actually make that second mesh valid; instead of an inverted element you'd have two overlapping elements! And I'm pretty sure you'll break some of our mesh topology assumptions (and thereby break the find_neighbors routine, leading to that remote_elem bug) in the process. --- Roy |
From: Joa L. <li...@jo...> - 2009-10-28 17:33:46
|
Hm... I almost giving up on this one... Looping through the elements, and moving nodes on the boundaries, it is not too difficult to make sure that the element I'm looking at isn't modified in a bad way. However, that node also belongs to a hole line of other elements, that might get destroyed. On top of that, when checking and correcting for hanging nodes, some of the elements that were good, just might go bad. I feel that if I move a node on the boundary given all the conditions I have to for fill I will not be able to improve very much over the "flat" surface approximation already in there. Maybe time to rethink the strategy: Since I produce the first mesh inside my code, using OpenCASCADE+Gmsh I'm thinking maybe I can let Gmsh regenerate a mesh for at each iteration, with the refined mesh of libmesh as a background mesh... or, I claim that I have regions that not are of interest... Why do I bother including them in the first place? cheers Joa On Mon, Oct 26, 2009 at 12:21:35PM -0500, Roy Stogner wrote: > > On Mon, 26 Oct 2009, Joa Ljungvall wrote: > > >What I'm trying to do is: > > > >1) Refine the mesh > >2) Find and store hanging nodes > >3) Move points on boundary to the "real" boundary > >then I do a do{...}While() including > >4) Check that Jacobian>0, if not switch node 0 and 2 (the pointers in the elem) > > I'm not sure this doesn't mess up something for neighbors.. > > This merely hides a symptom without fixing the problem. If you have > two elements ABC and CBD, and you invert one of them by moving a node > too far: > > A---------C A---------C > \ /| \ .-'/ > \ / | \ D / > \ / | \ | / > \ / | \|/ > B----D B > > Changing CBD into CDB doesn't actually make that second mesh valid; > instead of an inverted element you'd have two overlapping elements! > And I'm pretty sure you'll break some of our mesh topology assumptions > (and thereby break the find_neighbors routine, leading to that > remote_elem bug) in the process. > --- > Roy |
From: Roy S. <roy...@ic...> - 2009-10-28 17:38:50
|
On Wed, 28 Oct 2009, Joa Ljungvall wrote: > Since I produce the first mesh inside my code, using OpenCASCADE+Gmsh I'm > thinking maybe I can let Gmsh regenerate a mesh for at each iteration, with > the refined mesh of libmesh as a background mesh... If by "each iteration" you mean each refinement step in a steady solve, then this is likely to be the best solution for you. If you mean each time step in a transient solve, then regenerating the entire mesh means you'll now have projection error to worry about - how significant that is depends on your application and your formulation, though. --- Roy |
From: Roy S. <roy...@ic...> - 2009-10-14 22:15:37
|
On Wed, 14 Oct 2009, Derek Gaston wrote: > I thought I ran into similar problems when I did mesh redistribution / > smoothing + adaptivity. I was able to go in and "fix" the hanging > nodes somehow... But I don't remember where I did it ( I thought it > was down inside the library and I committed it back). Heh, and I thought your solution to that was to just do the adaptivity first and the redistribution second. ;-) > I'm not near a computer right now where I can go look at what I did > (I'm traveling this week)... But I bet if you look back a few years > (2007?) you might see a checkin from me referencing it. Aha! MeshTools::find_hanging_nodes_and_parents() looks like it's exactly the trick we'd need. Call that function, get the hanging_nodes map. Do hanging_node_moved = false Loop over hanging_nodes If a node isn't already in between its parents move it hanging_node_moved = true While(hanging_node_moved == true) The loop probably isn't the most efficient way to handle recursive constraints (maybe a depth-first descent into hanging_nodes at each entry?) but it would work. --- Roy |
From: Joa L. <li...@jo...> - 2009-10-15 07:45:07
|
Hi, Yes, it looks like the thing we need, expect that it only works for QUAD4... But I guess I could use it as a model and write a function in my program. If it works we add it to the svn version. Thanks for all the help. I think I understood the problem and how to fix it. I'll get back with news on the progress and/or more questions cheers Joa On Wed, Oct 14, 2009 at 05:15:26PM -0500, Roy Stogner wrote: > > On Wed, 14 Oct 2009, Derek Gaston wrote: > > >I thought I ran into similar problems when I did mesh redistribution / > >smoothing + adaptivity. I was able to go in and "fix" the hanging > >nodes somehow... But I don't remember where I did it ( I thought it > >was down inside the library and I committed it back). > > Heh, and I thought your solution to that was to just do the adaptivity > first and the redistribution second. ;-) > > >I'm not near a computer right now where I can go look at what I did > >(I'm traveling this week)... But I bet if you look back a few years > >(2007?) you might see a checkin from me referencing it. > > Aha! MeshTools::find_hanging_nodes_and_parents() looks like it's > exactly the trick we'd need. > > Call that function, get the hanging_nodes map. > > Do > hanging_node_moved = false > Loop over hanging_nodes > If a node isn't already in between its parents > move it > hanging_node_moved = true > While(hanging_node_moved == true) > > The loop probably isn't the most efficient way to handle recursive > constraints (maybe a depth-first descent into hanging_nodes at each > entry?) but it would work. > --- > Roy |
From: Joa L. <li...@jo...> - 2009-10-26 17:13:09
|
Hi, I've tried to get this running, but still have some problems... What I'm trying to do is: 1) Refine the mesh 2) Find and store hanging nodes 3) Move points on boundary to the "real" boundary then I do a do{...}While() including 4) Check that Jacobian>0, if not switch node 0 and 2 (the pointers in the elem) I'm not sure this doesn't mess up something for neighbors.. 5) Move the hanging nodes 6) If all Jacobian's>0 and all hanging nodes are in between there parents, go on... This works for a few iterations, then when looking for hanging nodes, i.e. the first thing I do after coursing and refining there is an element that has a neighbor with level==0 but that is not real... In the debugger Refining the mesh... [0] /usr/local/libmesh/include/geom/remote_elem.h, line 94, compiled Oct 11 2009 at 15:02:42 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Program received signal SIGABRT, Aborted. 0x00007fff80dcdff6 in __kill () (gdb) up #1 0x00007fff80e6f072 in abort () (gdb) #2 0x00007fff811185d2 in __gnu_cxx::__verbose_terminate_handler () (gdb) #3 0x00007fff81116ae1 in __cxxabiv1::__terminate () (gdb) #4 0x00007fff81116b16 in std::terminate () (gdb) #5 0x00007fff81116bfc in __cxa_throw () (gdb) #6 0x000000010282d7b4 in RemoteElem::n_sides () (gdb) #7 0x000000010189880e in Elem::which_neighbor_am_i (this=0x101903b80, e=<value temporarily unavailable, due to optimizations>) at elem.h:1293 1293 for (unsigned int s=0; s<this->n_neighbors(); s++) (gdb) print *neigh No symbol "neigh" in current context. (gdb) up #8 0x0000000101889418 in AGATAGeFEM::libMESH::find_hanging_nodes_and_parents (mesh=<value temporarily unavailable, due to optimizations>, hanging_nodes=@0x7fff5fbfd4f0) at ../../libs/libmesh/LaplaceProblem.cc:339 warning: Source file is more recent than executable. 339 neigh->which_neighbor_am_i(ancestor); (gdb) print *neigh $1 = { <ReferenceCountedObject<Elem>> = { <ReferenceCounter> = { _vptr$ReferenceCounter = 0x102d0ad90, static _n_objects = { _val = 605251 }, static _mutex = {<No data fields>} }, <No data fields>}, <DofObject> = { <ReferenceCountedObject<DofObject>> = { <ReferenceCounter> = { _vptr$ReferenceCounter = 0x102d0af08, static _n_objects = { _val = 605251 }, static _mutex = {<No data fields>} }, <No data fields>}, members of DofObject: old_dof_object = 0x0, static invalid_id = 4294967295, static invalid_processor_id = 65535, _id = 4294967295, _processor_id = 65535, _n_systems = 0 '\0', _n_v_comp = 0x0, _dof_ids = 0x0 }, members of Elem: static type_to_n_nodes_map = {2, 3, 4, 3, 6, 4, 8, 9, 4, 10, 8, 20, 27, 6, 15, 18, 5, 2, 4, 6, 8, 16, 18, 6, 16, 1, 0}, _nodes = 0x0, _neighbors = 0x0, _parent = 0x0, _children = 0x0, _rflag = 1 '\001', _pflag = 1 '\001', _p_level = 0 '\0', _sbd_id = 0 '\0', static _bp1 = 65449, static _bp2 = 48661 } So, any ideas? Bugs in the my implementation of logical issues? cheers Joa On Wed, Oct 14, 2009 at 05:15:26PM -0500, Roy Stogner wrote: > > On Wed, 14 Oct 2009, Derek Gaston wrote: > > >I thought I ran into similar problems when I did mesh redistribution / > >smoothing + adaptivity. I was able to go in and "fix" the hanging > >nodes somehow... But I don't remember where I did it ( I thought it > >was down inside the library and I committed it back). > > Heh, and I thought your solution to that was to just do the adaptivity > first and the redistribution second. ;-) > > >I'm not near a computer right now where I can go look at what I did > >(I'm traveling this week)... But I bet if you look back a few years > >(2007?) you might see a checkin from me referencing it. > > Aha! MeshTools::find_hanging_nodes_and_parents() looks like it's > exactly the trick we'd need. > > Call that function, get the hanging_nodes map. > > Do > hanging_node_moved = false > Loop over hanging_nodes > If a node isn't already in between its parents > move it > hanging_node_moved = true > While(hanging_node_moved == true) > > The loop probably isn't the most efficient way to handle recursive > constraints (maybe a depth-first descent into hanging_nodes at each > entry?) but it would work. > --- > Roy |
From: Roy S. <roy...@ic...> - 2009-10-14 13:59:37
|
On Wed, 14 Oct 2009, Joa Ljungvall wrote: > When I loop over the nodes, find a face that is on the boundary, and > then move the nodes the problem is, if I understood it correct, that > I might move a hanging node out of the line between the nodes of the > not so refined neighbour. This is not the only possible failure case: You also might move a vertex node without moving a hanging node (even a non-boundary hanging node) that depends on it. That's the only failure mode in 2D, in fact. > But I can check if this is the case, and instead of moving only the > hanging node I move all three nodes keeping them on a line to > minimize the error vis-a-vis the geometry. Or move all 6 nodes, in case you moved node A which depends on node B and C, node C depends on D and E, and node F depends on nodes B and G? In 3D this is not an extremely unlikely possibility, even with tets. That's why I suggested just fixing things with the constraint equations. It sounds great to optimize the geometry error, but one step at a time. > If this is a region where the solution changes fast this not so > refined element will be refined allowing an improved description of > the geometry. Also, in 3D "this not so refined element" may be "these not so refined elements", since any number of tets can share a boundary node or a boundary edge. --- Roy |
From: Joa L. <li...@jo...> - 2009-10-14 15:12:31
|
So returning to the "fix" using the constraint equations... I don't quite understand the idea behind this. If my hanging node no longer is at the right position, my mesh is broken, no? So fixing this with the constraint equations is a way of making the calculations go on, assuming that the error we introduced this way is small, and not a complete solution to the problem? This is way I came up with the idea of moving the nodes in groups... Joa Joa On Wed, Oct 14, 2009 at 08:38:59AM -0500, Roy Stogner wrote: > > On Wed, 14 Oct 2009, Joa Ljungvall wrote: > > >When I loop over the nodes, find a face that is on the boundary, and > >then move the nodes the problem is, if I understood it correct, that > >I might move a hanging node out of the line between the nodes of the > >not so refined neighbour. > > This is not the only possible failure case: You also might move a > vertex node without moving a hanging node (even a non-boundary hanging > node) that depends on it. That's the only failure mode in 2D, in > fact. > > >But I can check if this is the case, and instead of moving only the > >hanging node I move all three nodes keeping them on a line to > >minimize the error vis-a-vis the geometry. > > Or move all 6 nodes, in case you moved node A which depends on node B > and C, node C depends on D and E, and node F depends on nodes B and G? > In 3D this is not an extremely unlikely possibility, even with tets. > > That's why I suggested just fixing things with the constraint > equations. It sounds great to optimize the geometry error, but one > step at a time. > > >If this is a region where the solution changes fast this not so > >refined element will be refined allowing an improved description of > >the geometry. > > Also, in 3D "this not so refined element" may be "these not so refined > elements", since any number of tets can share a boundary node or a > boundary edge. > --- > Roy |