Thank you JP, your hints are highly appreciated. I will work on creating my own data structure, to see how complex and effective it can be. In the meanwhile though I made some tests with "almost" zero volumes cohesive interfaces and they are indeed simple and effective. So, if there is a chance to allow designing these elements in deal.ii I guess it is beneficial. I became aware that in many cases of interface problem using zero volume elements is a natural choice and papers have been extensively written on how to implement them in Abaqus or other commercial tools. If I can be of any help in promoting the discussion on github I will definitely do.
Thank you again, Alberto *Alberto Salvadori* Dipartimento di Ingegneria Civile, Architettura, Territorio, Ambiente e di Matematica (DICATAM) Universita` di Brescia, via Branze 43, 25123 Brescia Italy tel 030 3711239 fax 030 3711312 e-mail: alberto.salvad...@unibs.it web-pages: http://m4lab.unibs.it/faculty.html http://dicata.ing.unibs.it/salvadori On Mon, Aug 27, 2018 at 8:27 AM Jean-Paul Pelteret <jppelte...@gmail.com> wrote: > Hi Alberto, > > This is simply not possible. The deal.II triangulation contains as much > information as it needs, and we do not associate a normal with a volume. > Normals are computed internally (and not stored!) based not only on the > geometry of each real but also the mapping associated with the cell face. > So if we read this data in, and then you associated a non-flat manifold > with a face then you’d end up with a mis-match between the read-in normal > and that determined by the manifold definition and mapping. > > If it were possible to create a zero-thickness element, then one might > simply create a rule as to what you consider the “normal” to be. So, when > looping over the faces of the zero-thickness element (as one would do to > perform surface integration) then one might simply only perform the > integration on faces shared with, say, domain 1. This way the normal is > well defined, and completely independent of how the grid is created. > > Best, > J-P > > On 26 Aug 2018, at 12:07, Alberto Salvadori <alberto.salvad...@unibs.it> > wrote: > > As a short follow up of the former question, I wonder how to pass > geometrical information to a triangulation, if they are not included in the > file read by grid.in(). As an example connected to the problem described > above, consider the notion of "normal" at point 1/2 (in 3D: normal at the > surface that joins the two subdomains). One can easily identify a normal > for each "zero-thickness" element based on the gmsh file numeration, but it > looks to me that in the read.msh function the elm_number plays no role in > the triangulation defined by deal.ii. If this is correct, how to pass the > geometrical info to the deal.ii triangulation? > Thank you very much for any suggestions, > Alberto > > Il giorno venerdì 24 agosto 2018 14:39:41 UTC+2, Alberto Salvadori ha > scritto: >> >> Thank you Wolgang and Jean Paul. >> >> In these days I made some implementation in the line you propose but I am >> not happy with it and thus I am asking advices on a different path of >> reasoning. >> >> 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?). >> >> 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. >> >> I appreciate your comments. >> Best >> >> Alberto >> >> >> >> 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> >>> 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> >>> 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> >>> (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 >>> www: >>> http://www.math.colostate.edu/~bangerth/ >>> >> > > Informativa sulla Privacy: http://www.unibs.it/node/8155 > > -- > 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. > > > -- > 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. > -- 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/ 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.
