projecting values of the solution on the boundary

[Alberto Salvadori]

Dear community,

I am asking some advices on the following issue. I am solving a simple problem, say a Poisson problem in the unknown field u, and a more involved problem separately. This second problem requires the values of u in the Neumann boundary conditions.

Accordingly, I guess one could solve the Laplacian first and calculate the numerical solution for u, say un. Afterwards one builds a solver for the more complex operator and in the Neumann part of the code - that may look like this for parallel::shared triangulations:

  for (unsigned int face_number=0; face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)

        if (

            cell->face(face_number)->at_boundary()

            &&

            cell->face(face_number)->boundary_id() == 2    // Neumann boundaries

            )

        {

          

          fe_face_values.reinit (cell, face_number);

          

          // define points and normals

          

          std::vector< Point<dim> >    points  = fe_face_values.get_quadrature_points();

          std::vector< Tensor<1,dim> > normals = fe_face_values.get_all_normal_vectors();

          

          // calculate neumann values

          

          for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)

          {

            

            // values: mechanical

            Tensor<1,dim> mech_neumann_value;

            neumann_bc_for_mech.bc_value( points[q_point], normals[q_point], mech_neumann_value );

             .....

     

one needs the value of the field un at points  points[q_point] to be passed to neumann_bc_for_mech.bc_value.

Which is an effective way to calculate this amount ?

Thank you.

Alberto

[Jean-Paul Pelteret]

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() 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() function so you could create a lookup between a point in problem 1 and a cell in problem 2 again using GridTools::Cache (although I’m guessing that this is what it does internally).

I hope that this gives you some ideas to play with!

Best,

J-P

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+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

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/

[Alberto Salvadori]

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

[Alberto Salvadori]

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

[Jean-Paul Pelteret]

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

[Jean-Paul Pelteret]

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

[Alberto Salvadori]

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....@unibs.it
web-pages:
 http://m4lab.unibs.it/faculty.html
 http://dicata.ing.unibs.it/salvadori