Hi Alberto, > In brief: I'd like to design an algorithm to solve a problem over two > separate domains connected by an interface. Think in 1D of line 0_1 > separated in two parts 0_1/2 and 1/2_1. > In the simplest case, a linear elastic problem on each subdomain with a > spring that joins the two sides located at point 1/2. The boundary conditions > of the two parts at point 1/2 are a linear combination of the solutions on > the two boundaries. > The weak form of this problem can easily be written and surface integrals > involving unknowns arise. I tried to solve this with the algorithm suggested > by Wolfgang, and it works fine for very small meshes. However, this strategy > would imply that for each cell on the boundary (at point 1/2 for the domain > 0_1/2) one has to seek trough the whole triangulation for the cell that > corresponds to point 1/2 in the remaining part (1/2_1). This procedure is > very expensive, or at least that comes out to me. (By the way, although the > triangulation is parallel shared I can only see cells on the same node: shall > I use a specific iterator?).
Have you considered creating your own data structure that represents the cohesive zone, and precomputing these inter-domain relations? What I would consider doing is something along the lines of what you’d do in contact mechanics. You could maybe mark one domain boundary as the “master” surface and the other the “slave”, and then find out and define once up front (or each time you refine/coarsen) the relationship between the faces of the two domains in the cohesive zone. You could then do all of standard volume and face assembly independent of the cohesive zone, and then loop over all faces in the cohesive zone (as stored in your data structure, not via the triangulation / dofhandler) and do the final integration on this subdomain. In this way you don’t recompute the expensive geometric problem each time you do assembly. > To circumvent the issue, one could define a zero-thickness element at point > 1/2 (called cohesive in the literature sometimes) because in such a case the > element brings in all the connectivities that are required and one can still > define the tangent stiffness matrix based on surface integrals. I have been > working on this idea but I run into a major problem in loading the > triangulation from gmsh, since elements with zero volume are not considered > to be good as for now. I wonder if one could get rid of this control or if > zero-volume condition is used all over deal.ii as a check of something going > wrong. > > In fact one could also move the nodes a bit to generate "almost" zero > volumes. In some cases it can be done easily, but in general this approach is > unfeasible for complex meshes. Yeah, this is not going to work. There’s a check performed when reading in a triangulation that all cells have a positive volume. This makes sense in the overwhelming majority of use cases, so I’m quite confident that we wouldn’t consider removing it. Maybe we’d consider having a flag sent in to the read_mesh() functions that disables the check (or some equivalent mechanism), but I think thats a discussion worth having on github rather than here. Best regards, J-P > > Il giorno lunedì 13 agosto 2018 21:47:29 UTC+2, Wolfgang Bangerth ha scritto: > On 08/08/2018 02:03 AM, Jean-Paul Pelteret wrote: > > Hi Alberto, > > > > I have dealt with a similar problem where I needed to transmit solution > > values between two different problems that shared a common interface. If > > I remember correctly, the way that I did this was to precompute the > > positions at which I would need the solutions from problem 1 for problem > > 2. In a post processing step for problem 1 I then computed this data up > > front (cache it) and then fetch it from problem 2. This avoided the > > issue of going to look for which cell in problem 1 a point in problem 2 > > lay. I did this all manually, but I guess that this approach could be > > made simpler now that we have GridTools::compute_point_locations() > > <https://www.dealii.org/developer/doxygen/deal.II/namespaceGridTools.html#a8e8bb9211264d2106758ac4d7184117e > > > > <https://www.dealii.org/developer/doxygen/deal.II/namespaceGridTools.html#a8e8bb9211264d2106758ac4d7184117e>> > > which > > allows a speedy lookup between a point and its containing cell. > > > > FEFieldFunction has a cache so you could also investigate using it > > (maybe it wouldn’t be too slow for your case). It also has the > > set_active_cell() > > <https://www.dealii.org/9.0.0/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html#a0206a45c90d523792eea8bd725d14788 > > > > <https://www.dealii.org/9.0.0/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html#a0206a45c90d523792eea8bd725d14788>> > > function > > so you could create a lookup between a point in problem 1 and a cell in > > problem 2 again using GridTools::Cache > > <https://www.dealii.org/developer/doxygen/deal.II/classGridTools_1_1Cache.html > > > > <https://www.dealii.org/developer/doxygen/deal.II/classGridTools_1_1Cache.html>> > > (although > > I’m guessing that this is what it does internally). > > This is indeed what you would do if you had separate meshes. The better > approach, of course, is to use the same mesh for both problems. In that > case, you can use this to obtain the gradients of the Poisson solution > un at the quadrature points on the faces where you need these values for > the second problem: > FEFaceValues poisson_fe_face_values (...); > std::vector<Tensor<1,dim>> un_gradients (...); > std::vector<double> un_dot_n_values (...); > > for (cell=...) > for (face=...) > if (face is interesting) > { > poisson_fe_face_values.reinit (cell, face); > poisson_fe_face_values.get_function_gradients (poisson_solution, > un_gradients); > for (q=0...) > un_dot_n_values[q] = un_gradients[q] * > poisson_fe_face_values.normal_vectors(q); > > ...do assembly for the second problem using the values of > (grad un).n just computed... > } > > Best > W. > > > -- > ------------------------------------------------------------------------ > Wolfgang Bangerth email: bang...@ <>colostate.edu > <http://colostate.edu/> > www: http://www.math.colostate.edu/~bangerth/ > <http://www.math.colostate.edu/~bangerth/> > > > Informativa sulla Privacy: http://www.unibs.it/node/8155 > <http://www.unibs.it/node/8155> > > -- > 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>. > For more options, visit https://groups.google.com/d/optout > <https://groups.google.com/d/optout>. -- 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. For more options, visit https://groups.google.com/d/optout.
