HI Manuel, I’ve got a brief answer to your followup questions, but let me first say that Andrew’s suggestion, which seems to be the most elegant/optimal approach to the problem, is one that is possible to for you to realise. It’s also made me think that perhaps some of the tools that have been added somewhat recently to help tackle immersed problems could be of help? I have zero experience with these tools, so I’m just going to put it out there for you to investigate and draw your own conclusions from:
https://dealii.org/current/doxygen/deal.II/namespaceNonMatching.html <https://dealii.org/current/doxygen/deal.II/namespaceNonMatching.html> https://dealii.org/current/doxygen/deal.II/step_60.html <https://dealii.org/current/doxygen/deal.II/step_60.html> I’m sure that if you have any questions on those, then the other developers would be willing to help out to answer them! > My question is what happens to the DoFs that are shared by a cell with a > FE_Nothing space and another cell with a FE_Q space? This is of course going > to happen at all cells in the boundary since some DoFs of those cells are > indeed in the boundary while others are still internal DoFs. Does this make > sense? I suppose I can apply a constraint to those internal DoFs, similar to > my previous "hack”. Yes, that’s correct. You’d be doing exactly the same thing as you’d detailed before, but just for fewer elements / DoFs. > The advantage following your suggestion is that the number of DoFs that I > need to constrain is much smaller. Indeed. So not only is your global linear system smaller, but the AffineConstraints object also has less entries. Why this might be important is due to a comment that can be found in the AffineConstraints introductory notes <https://dealii.org/current/doxygen/deal.II/classAffineConstraints.html>: The class is meant to deal with a limited number of constraints relative to the total number of degrees of freedom, for example a few per cent up to maybe 30 per cent; and with a linear combination of M other degrees of freedom where M is also relatively small (no larger than at most around the average number of entries per row of a linear system). It is not meant to describe full rank linear systems. Its hard to say whether this might have a performance effect wrt. what you’re trying to do, but nevertheless using Fe_Nothing to limit the number of DoFs that you have to constrain would likely mean that you could consider finer discretisation before you have to start considering what the impact of the number of constraints have on performance (and anything else). Best, Jean-Paul > On 29. Mar 2021, at 22:27, manuel.q...@gmail.com > <manuel.quezada...@gmail.com> wrote: > > Hi Jean-Paul, > > I have taken a look at FE_Nothing and hp::DoFHandler. If I understand > correctly, the idea is that all the interior cells will use an FE_Nothing > space and the boundary cells will use a standard FE_Q space, right? My > question is what happens to the DoFs that are shared by a cell with a > FE_Nothing space and another cell with a FE_Q space? This is of course going > to happen at all cells in the boundary since some DoFs of those cells are > indeed in the boundary while others are still internal DoFs. Does this make > sense? I suppose I can apply a constraint to those internal DoFs, similar to > my previous "hack". The advantage following your suggestion is that the > number of DoFs that I need to constrain is much smaller. > > Thanks for your suggestion. I think this is a very nice approach, > > Manuel > > > On Monday, March 29, 2021 at 8:09:39 PM UTC+3 Jean-Paul Pelteret wrote: > Hi Manuel, > > As a small addition to what Wolfgang has already suggested, you can also > considering using the hp framwwork in conjunction with FE_Nothing > <https://dealii.org/current/doxygen/deal.II/classFE__Nothing.html> elements > in order to, effectively, disable elements in the interior of the volumetric > domain instead of constraining them directly. This is the equivalent of > saying that you simply wouldn’t assign DoFs for those elements, but rather > only for the elements that are positioned along the boundary of the domain. > You’d still have to constrain the DoFs that are off of the surface itself, > but at least there are a lot fewer DoFs that require this. > > Best, > Jean-Paul > > >> On 29. Mar 2021, at 18:59, Wolfgang Bangerth <bang...@colostate.edu >> <applewebdata://82C121CE-670B-46E8-827C-F024382CF28F>> wrote: >> > >> >> Manuel, >> >>> I need to solve a system of many PDEs (currently 18). The actual number >>> changes since the PDEs are defined by a truncated infinite sum and some >>> recursive equations. Anyway, the first two equations are defined on a 2D or >>> 3D domain (call it Omega) and the rest is only defined on the boundary of >>> Omega. My first idea was to create a FESystem object and attach to it 2 FE >>> spaces defined on Omega and 16 FE spaces defined on the boundary. That is, >>> I tried these options >>> * FESystem<1, 2> fe(FE_Q<2, 2>(degree), 2, FE_Q<1, 2>(degree), 16) >>> * FESystem<2, 2> fe(FE_Q<2, 2>(degree), 2, FE_Q<1, 2>(degree), 16) >>> However, that produced compilation errors. What I understand is that the >>> FESystem must be created by either FE spaces defined in Omega or FE spaces >>> defined in the manifold, but not a combination. >> >> Yes, that is correct -- because FESystem<dim,spacedim> is derived from >> FiniteElement<dim,spacedim>, the elements that form the system must have a >> unique 'dim' parameter. >> >> >>> As a way around I tried (a very ugly hack) where I defined all FE spaces in >>> Omega and then add constraints to all DoFs associated with the internal >>> nodes for the variables that are defined on the manifold. This sounds like >>> a tremendous waste of resources since most of the equations in my system >>> live in the manifold, so I end up constraining most of the DoFs. And even >>> worse, it produced a matrix that seems to be sometimes singular and >>> sometimes non-singular. That is, I can solve the system and get accurate >>> results (I have the exact solution) for some grids but in other grids the >>> matrix is singular so I get NaNs. I suspect this problem is a consequence >>> of constraining most of my DoFs. FYI, the solver that I am using is the >>> MUMPS implementation in PETSc. >> >> I don't know about the origins of the illposedness, but this is a valid >> approach. ASPECT, for example, does this kind of thing when implementing a >> free boundary. You could take a look here >> https://github.com/geodynamics/aspect/blob/master/include/aspect/mesh_deformation/free_surface.h >> >> <https://github.com/geodynamics/aspect/blob/master/include/aspect/mesh_deformation/free_surface.h> >> https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc >> >> <https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc> >> >> In particular, ASPECT just doesn't do anything at all with the interior >> degrees of freedom in the solve phase: >> https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc#L345-L352 >> >> <https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc#L345-L352> >> In other words, we happily solve a linear system with a matrix with a lot of >> zero rows and columns. I've got questions about this here: >> https://github.com/geodynamics/aspect/issues/3110 >> <https://github.com/geodynamics/aspect/issues/3110> >> But it seems to work. >> >> >>> I wonder if there is a proper or standard way to construct systems that are >>> defined by some equations in a domain and some in a manifold of the domain. >>> I definitely don't think constraining most of my DoFs is the way to go. >> >> The "proper" way to deal with this would probably be to have two DoFHandler >> objects with dim and dim-1, and then to build a coupled linear system. I >> believe that we have most of the building blocks for this already in place >> (take, for example, a look at step-60 and step-70) but there is no concrete >> example of coupling bulk and surface problems. It has been on my TODO list >> for quite a while already to write a tutorial program for this and see what >> it would take to make that work properly. >> >> Out of curiosity, what is the problem you are trying to solve? >> >> Best >> W> >> >> -- >> ------------------------------------------------------------------------ >> Wolfgang Bangerth email: bang...@colostate.edu >> <applewebdata://82C121CE-670B-46E8-827C-F024382CF28F> >> www: http://www.math.colostate.edu/~bangerth/ >> <http://www.math.colostate.edu/~bangerth/> >> > >> -- >> The deal.II project is located at http://www.dealii.org/ >> <http://www.dealii.org/> >> For mailing list/forum options, see >> https://groups.google.com/d/forum/dealii?hl=en >> <https://groups.google.com/d/forum/dealii?hl=en> >> --- You received this message because you are subscribed to the Google >> Groups "deal.II User Group" group. >> To unsubscribe from this group and stop receiving emails from it, send an >> email to dealii+un...@googlegroups.com >> <applewebdata://82C121CE-670B-46E8-827C-F024382CF28F>. >> To view this discussion on the web visit >> https://groups.google.com/d/msgid/dealii/6bb61a54-9428-2174-48ae-83b649c1e3a7%40colostate.edu >> >> <https://groups.google.com/d/msgid/dealii/6bb61a54-9428-2174-48ae-83b649c1e3a7%40colostate.edu>. > > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <https://groups.google.com/d/forum/dealii?hl=en> > --- > You received this message because you are subscribed to the Google Groups > "deal.II User Group" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to dealii+unsubscr...@googlegroups.com > <mailto:dealii+unsubscr...@googlegroups.com>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/b09a005d-8849-4c70-a851-1997aa398ec5n%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/b09a005d-8849-4c70-a851-1997aa398ec5n%40googlegroups.com?utm_medium=email&utm_source=footer>. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/82CB463B-B8F1-4F3A-879C-1B5910ACBA62%40gmail.com.
