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 <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.
