Re: [deal.II] projecting values of the solution on the boundary

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

Re: [deal.II] projecting values of the solution on the boundary

[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

> On 26 Aug 2018, at 12:07, Alberto Salvadori  
> 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() 
> >  >  
> > 

Re: [deal.II] projecting values of the solution on the boundary

[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

> 
> 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() 
> >  >  
> > >
> >  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 
> > 

Re: [deal.II] projecting values of the solution on the boundary

[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

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 

Re: [deal.II] projecting values of the solution on the boundary

[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



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> un_gradients (...); 
>std::vector 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: 

Re: [deal.II] projecting values of the solution on the boundary

[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() 
 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).


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> un_gradients (...);
  std::vector 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: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

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

Re: [deal.II] projecting values of the solution on the boundary

[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

> On 07 Aug 2018, at 17:47, Alberto Salvadori  
> wrote:
> 
> 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::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 >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 
>   {
> 
> 
> // 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
> 
> 
> 
> 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.

[deal.II] 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::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 >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 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


-- 


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.