from:"Sylvain Mathonnière"

Re: [deal.II] Loss of symmetry by changing parameter

2021-10-12 Thread Sylvain Mathonnière

For the record, the lines of code in my last post were wrong. It should be 
for the boundary worker :

// new cell iterator for other dof_handler 
typename DoFHandler::active_cell_iterator cell_2 (>
get_triangulation(), cell->level(), cell->index(), _handler_HE);

// Initialise FEValues with FACE quadrature formula and reinitialise 
with current cell 
FEFaceValues fe_values_HE_bound(mappingQ, fe_HE, 
quadrature_face_formula, update_values);
fe_values_HE_bound.reinit(cell_2, face_no);

std::vector T_old_data_bound(q_points.size());
fe_values_HE_bound.get_function_values(old_T,T_old_data_bound);

Best,

Sylvain M.

El viernes, 8 de octubre de 2021 a la(s) 14:30:54 UTC+2, Wolfgang Bangerth 
escribió:

> On 10/8/21 5:13 AM, Sylvain Mathonnière wrote:
> > How can one suspect there is such problem with the code in those cases ? 
> Just 
> > with experience I guess and doing multiple simulations ?
>
> Making every mistake myself at least once in the last 20 years has 
> certainly 
> helped me have an eye for things :-)
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1fcdc32c-3943-4c49-b8bc-f5a218f0e2f0n%40googlegroups.com.

Re: [deal.II] Loss of symmetry by changing parameter

2021-10-08 Thread Sylvain Mathonnière

Oh thank you very much ! That was indeed the problem, I changed the code by 
calling FEValues using the face quadrature formula and used it locally to 
avoid those race conditions. It looks like this :

// new cell iterator for other dof_handler
typename DoFHandler::active_cell_iterator cell_2 (>
get_triangulation(), cell->level(), cell->index(), _handler_HE);

// Initialise FEValues with FACE quadrature formula and reinitialise 
with current cell
FEValues fe_values_HE(mappingQ, fe_HE, q_points, update_values);
fe_values_HE.reinit(cell_2);

// get data at quadrature point ON THE FACE
std::vector T_old_data(q_points.size());
fe_values_HE.get_function_values(old_T,T_old_data);

You managed to solve 2 of the problem I was having. I believe the code 
works as intended now.
Those kind of errors are quite tricky to catch. It did not show much in the 
results I was getting and I am pretty sure the solution I had (except the 
set parameter 3 which was obviously wrong) would have converged by 
decreasing the size of the mesh.
Moreover the method of manufactured solution could not catch this error at 
all since it was in the RHS.
How can one suspect there is such problem with the code in those cases ? 
Just with experience I guess and doing multiple simulations ?

Regardless, thank you very much again ! I hope one day I will be able to 
help back this community as much as it helped me so far.

Best,

Sylvain Mathonnière

El jueves, 7 de octubre de 2021 a la(s) 23:04:15 UTC+2, Wolfgang Bangerth 
escribió:

> On 10/7/21 8:36 AM, Sylvain Mathonnière wrote:
> > 
> > The main issue in the code is the assembly, more precisely on line 305 : 
> > copy_data.cell_rhs_RTE(i) +=  pow(T_old_data[point],4) * 1e-8 * 
> > fe_face.shape_value(i, point)*JxW[point];
> > if I remove the pow(T_old_data[point],4), then it gets perfectly 
> > symmetric but otherwise no. I cannot get my mind around the reason why.
>
> Let's do some more remote debugging then :-) Since you've already 
> identified that the problem appears to be the old values T_old_data, 
> output those as well. You should get ... actually, I found your bug: You 
> are looping over the quadrature points of the FEFaceValues object here:
>
> const auto _points = fe_face.get_quadrature_points();
> ...
> for (unsigned int point = 0; point < q_points.size(); ++point)
> {
> double beta_dot_n = direction_temp * normals[point];
>
> if (beta_dot_n < 0)
> {
> for (unsigned int i = 0; i < n_facet_dofs; ++i)
> {
> // problem is here
> copy_data.cell_rhs_RTE(i) += pow(T_old_data[point],4) * 1e-8
> * fe_face.shape_value(i, point)*JxW[point];
>
> }
>
> So one would expect that T_old_data[point] is also indexed by the 
> quadrature points on the face. But you initialize
>
> fe_values_HE.reinit(cell_2);
> fe_values_HE.get_function_values(old_T,T_old_data);
>
> that is, T_old_data is indexed by the quadrature points on the *cell*, 
> not the face. This can't be right.
>
> Separately, you use the same T_old_data object for all face_workers. But 
> these may be called in parallel, and you will get race conditions. You 
> need to make sure you have a separate T_old_data object on every thread, 
> for example by making the variable local to the face_worker lambda 
> function.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/b526b942-7ef7-42ca-b6d5-d8ae25849d94n%40googlegroups.com.

[deal.II] Re: Loss of symmetry by changing parameter

2021-10-06 Thread Sylvain Mathonnière

Always some errors are left after I wrote the message... errare humanum est.
I inversed the figure in the text, I meant : 

"On the fig. 5, the bottom right cell has its form function different than 
0 for dof =0 (angle dof) (i.e *fe_face.shape_value(dof = 0; point_gauss) 
=/= 0* ) 
whereas on the fig. 4, the bottom left cell has *fe_face.shape_value(dof = 
0; point_gauss) *= 0 for all gauss point.

Also when I wrote "The BC", I wanted to say that the boundary conditions 
are the same in both cases and are actually implemented directly in the 
assembly (like in step-12).

El miércoles, 6 de octubre de 2021 a la(s) 17:36:29 UTC+2, Sylvain 
Mathonnière escribió:

> Dear all,
>
> I am experiencing a symmetry issue in my calculation with certain set of 
> parameters.
> In the 3 figures below, I used different sets of parameters (same 
> equation) and plotted the result once the simulation reaches steady state.
> The two first figures look fairly symmetric (as it should be) but the 3rd 
> one is obviously not symmetric.
> Upon closer inspection, one can see that the bottom left corner looks a 
> bit "shaky" on the first 2 figures.
>
> [image: fig1_3.PNG]
>
> I tried to get to the bottom of why this happens numerically ? 
> I found that at some point during my assembly of the RHS, there is a 
> difference in what I obtained even though it should be perfectly symmetric. 
> The solver takes also significantly longer to solve it.
> I represented the RHS on a coarse grid to highlight the difference.
>
> [image: fig4_5.PNG]
>
> By digging even more I realised that the discrepancy steams from the 
> boundary_worker (similar to the one of step-12).
> On the fig. 4, the bottom right cell has its form function different than 
> 0 for dof =0 (angle dof) (i.e *fe_face.shape_value(dof = 0; point_gauss) 
> =/= 0* ) 
> whereas on the fig. 5, the bottom left cell has *fe_face.shape_value(dof 
> = 0; point_gauss) *= 0 for all gauss point.
>
> For the record, I am solving several time the radiative transport equation 
> I  using the discontinuous galerkin method using several directions. 
> [image: RTE.PNG]
> So the fig.4 is the RHS for one direction beta=(0.14; 0.14) and fig 5. is 
> the RHS for the other direction beta=(0.14;-0.14). The BC 
>
> *My question is the following: *
> Is this a normal behaviour that I am experiencing and the lack of symmetry 
> in fig. 3 is linked with the fact that my mesh is not fine enough ?
> or could it be the result of something wrong in my code (I attached a 
> snipset of it below) ?
>
> I am not experienced enough to know the difference. I know the code works 
> well because I compared it to scientific reference + using the method of 
> manufactured solution, so if there is a mistake, it has to be "small".
>
> Best regards,
>
> Sylvain Mathonnière 

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/dda13b58-6234-45bd-8a29-3fbf0d7e6ae3n%40googlegroups.com.

Re: [deal.II] Re: function face_component_to_system_index() for Discontinuous Galerkin ?

2021-09-30 Thread Sylvain Mathonnière

 *About writing the function :*
I must confess I did not try to write the function myself but I did check 
to see how it was implemented and it seems like it is just retrieving the 
value from a table. I currently only have the user version of deal.II (v 
9.2.0),
If I wanted to write this function, the proper way to do it would be to 
download the developer version 9.3.0, write the function down then 
recompile the librairies ?
I haven't thought much about contributing before to deal.II because I am 
fairly new to it (I love it) and not that skilled at c++.

*About your first method :*
Yes, in the end this is the way I implemented it, but it takes really long 
time to compute since *dof_per_cell = 160* for me (with Q1 mapping and 
fe_degree =1), so I have 160*160 iterations but the "if condition" is only 
true 1/40th of the time, so I end up with a lot of unnecessary looping.
My equation system is *b.\nabla I + I = RHS* where *b* is a direction 
vector and I is a scalar function. For the discrete ordinate method, I must 
solve this equation for 40 values of *b*.

*About your second method :*
The problem in doing so, is that for me the 40 components are already 
included in the *dof_per_cell* (4 nodes * 40 components = 160 
dof_per_cell), so the ith and jth indices of the loops already correspond 
to a specific component. Hence why I used the  
system_to_component_index() function.
( I have  *FESystem * =  *fe_eq1 (FE_DGQ(1), 40) )* 

*What worked the best :*
I mentioned in my previous post that since my equations are not coupled 
with each other, I could technically loop over the assembly and solving 
functions 40 times only changing a constant every time (I change the 
direction vector  for the discrete ordinate method).
I though it was a really unefficient way but I tried it nonetheless and it 
turns out that it is actually much faster (I gain a speed factor of 16). 
In the end, assembling and solving 40 times a system with a sparse matrix 
of size (n_cells*4)^2 is faster than assembling and solving 1 time a system 
with a sparse matrix of size (n_cells*4*40)^2.
It sounds correct but does that make sense ?

*Caveat : *
I use MeshWorker::mesh_loop with queue_length = 1 so I believe I only use 
one thread, maybe it is a different result with more thread.
Otherwise if I use mesh_loop like in step 12, I get an error 
"pure virtual method called
terminate called without an active exception
Abandon (core dumped)"
I believe this error occurs because I create a second iterator in the cell 
worker to access values bound to another dof_handler like this: 
*typename DoFHandler::active_cell_iterator cell_2 
(>get_triangulation(), cell->level(), cell->index(), _handler_2);*

Best,

Sylvain

El miércoles, 29 de septiembre de 2021 a la(s) 20:37:05 UTC+2, Wolfgang 
Bangerth escribió:

> On 9/23/21 5:48 AM, Sylvain Mathonnière wrote:
> > 
> > I found in the documentation a fonction 
> > *face_system_to_component_index()* but no 
> > function*face_component_to_system_index()*. What would be the best way 
> > to proceed at this point ?
>
> It would not be difficult to actually add this function. Would you like 
> to give it a try? If you see how the non-face functions are implemented, 
> you'd probably see right away how to add such a function for faces as well.
>
> That said, the way you approach this is a bit unusual. A better way 
> would be to consider that you have a vector-valued element and, as we 
> always do, loop over all DoFs i and j that belong to *all* components. 
> You can then either do something like
>
> for (i=0...dofs_per_cell)
> for (j=0...dofs_per_cell)
> if (fe.system_to_component_index(i).first ==
> fe.system_to_component_index(j).first)
> ...add matrix term...
>
> or, if you would rather do this, write things as
>
> for (i=0...dofs_per_cell)
> for (j=0...dofs_per_cell)
> for (c=0...n_components)
> copy_data_face.cell_matrix_RTE(i,j) +=
> ... fe_face[c].shape_value(i,q) ...
> ... fe_face[c].shape_value(j,q) ...
>
> The first of these ways would correspond to how step-8 assembles its 
> matrix, whereas the second is used in all other tutorial programs 
> solving vector-valued problems (e.g., step-22). You might also want to 
> look at the vector-valued documentation module for discussions around 
> these sorts of issues.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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 em

[deal.II] Re: function face_component_to_system_index() for Discontinuous Galerkin ?

2021-09-23 Thread Sylvain Mathonnière

I wanted to add some clarification:

I realised that what I called *n_vertices* is missleading, it is not the 
number of vertices I want but the number of "nodes per cell". Numerically, 
it is however correct if  *FE_DGQ(1)*.
More generally, it should be something like *n_nodes = 
(fe_eq1.degree+1)^dim*, so n_nodes = 4 for  *FE_DGQ(1)* but n_nodes = 
9 for *FE_DGQ(2).*

El martes, 21 de septiembre de 2021 a la(s) 10:00:32 UTC+2, Sylvain 
Mathonnière escribió:

> Hello,
>
> I have an equations (named eq1) that I am solving just like in step-12 
> using the discontinuous Galerkin method but I am solving it using the 
> discrete ordinate method so I actually have *N_dir = 40* directions so 40 
> equations. 
> My setup looks like this :
>
> * FESystemfe_eq1; *
> * DoFHandler  dof_handler_eq1;*
>
> And the constructor of the class :
>
> * template  DG_FEM::DG_FEM ()*
> * : fe_eq1 (FE_DGQ(1),N_dir), dof_handler_eq1 (triangulation)) {}*
>
> For the eq1, I am using like in *step-12* the *MeshWorker::mesh_loop(...)* 
> to assemble my matrix :
>
>
> In the cell_worker function, in order to access the index_dof_i-th 
> component (out of 40) of my shape_function for the direction dir_comp, I 
> have:
>
> *...*
> *for (unsigned int dir_comp = 0; dir_comp < N_dir; ++dir_comp)*
> *{*
> * for (unsigned int i = 0; i < n_vertices; ++i) // here I hardcoded 
> n_vertices = 4 for quadrangle*
> * {*
> * unsigned int index_dof_i = fe_eq1.component_to_system_index(dir_comp,i);*
>
> * for (unsigned int j = 0; j < n_vertices; ++j)*
> * {*
> * unsigned int index_dof_j = fe_eq2.component_to_system_index(dir_comp,j);*
>   
> * copy_data_face.cell_matrix_RTE(index_dof_j, index_dof_i) +=*
> * // some calculations involving*
> * fe_face.shape_value_component(index_dof_i, point, dir_comp)*
> * fe_face.shape_value_component(index_dof_j, point, dir_comp)*
> * }*
> * }*
> *}*
> *...*
>
> I need to loop on the vertices of my cell to do the assembly and I did not 
> find a better way than that. I cannot do like in step-8 looping over all 
> the dof_per_cell, because then
> I could not couple the component at one vertice with the same component at 
> another vertice, like I am doing above.
> I guess it is fine as long as I only have quadrangle. 
>
>
>
> This code works fine it seems and I can apply the same strategy for the 
> boundary_worker. My main issue is that this technique does not work for the 
> face_worker. Here I need to loop on all the vertices twice (because there 
> are two cells sharing the same vertices so two shape function per vertices).
> If I write :
>
> for (unsigned int j = 0; j < n_vertices*2; ++j)
> {
> unsigned int index_dof_jfe_eq2.component_to_system_index(dir_comp,j);
> ...
> this will trigger an error since *component_to_system_index(dir_comp,j)* 
> will not accept j >= 4
>
> I found in the documentation a fonction *face_system_to_component_index()* 
> but no function* face_component_to_system_index()*. What would be the 
> best way to proceed at this point ?
>  
> Best,
>
> Sylvain M.
>
> PS : As I write this email I realise that I could technically solve 40 
> times the equation, only changing the direction each time since the 
> components are not coupled to each other (they are howerver coupled with 
> another equation not mentionned here) but that looks cumbersome at best and 
> probably inefficiency.
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2334bf8a-a13a-4ec2-875f-333974dd6f59n%40googlegroups.com.

[deal.II] function face_component_to_system_index() for Discontinuous Galerkin ?

2021-09-21 Thread Sylvain Mathonnière

Hello,

I have an equations (named eq1) that I am solving just like in step-12 
using the discontinuous Galerkin method but I am solving it using the 
discrete ordinate method so I actually have *N_dir = 40* directions so 40 
equations. 
My setup looks like this :

* FESystemfe_eq1; *
* DoFHandler  dof_handler_eq1;*

And the constructor of the class :

* template  DG_FEM::DG_FEM ()*
* : fe_eq1 (FE_DGQ(1),N_dir), dof_handler_eq1 (triangulation)) {}*

For the eq1, I am using like in *step-12* the *MeshWorker::mesh_loop(...)* 
to assemble my matrix :


In the cell_worker function, in order to access the index_dof_i-th 
component (out of 40) of my shape_function for the direction dir_comp, I 
have:

*...*
*for (unsigned int dir_comp = 0; dir_comp < N_dir; ++dir_comp)*
*{*
* for (unsigned int i = 0; i < n_vertices; ++i) // here I hardcoded 
n_vertices = 4 for quadrangle*
* {*
* unsigned int index_dof_i = fe_eq1.component_to_system_index(dir_comp,i);*

* for (unsigned int j = 0; j < n_vertices; ++j)*
* {*
* unsigned int index_dof_j = fe_eq2.component_to_system_index(dir_comp,j);*
  
* copy_data_face.cell_matrix_RTE(index_dof_j, index_dof_i) +=*
* // some calculations involving*
* fe_face.shape_value_component(index_dof_i, point, dir_comp)*
* fe_face.shape_value_component(index_dof_j, point, dir_comp)*
* }*
* }*
*}*
*...*
   
I need to loop on the vertices of my cell to do the assembly and I did not 
find a better way than that. I cannot do like in step-8 looping over all 
the dof_per_cell, because then
I could not couple the component at one vertice with the same component at 
another vertice, like I am doing above.
I guess it is fine as long as I only have quadrangle. 



This code works fine it seems and I can apply the same strategy for the 
boundary_worker. My main issue is that this technique does not work for the 
face_worker. Here I need to loop on all the vertices twice (because there 
are two cells sharing the same vertices so two shape function per vertices).
If I write :

for (unsigned int j = 0; j < n_vertices*2; ++j)
{
unsigned int index_dof_jfe_eq2.component_to_system_index(dir_comp,j);
...
this will trigger an error since *component_to_system_index(dir_comp,j)* 
will not accept j >= 4

I found in the documentation a fonction *face_system_to_component_index()* 
but no function* face_component_to_system_index()*. What would be the best 
way to proceed at this point ?
 
Best,

Sylvain M.

PS : As I write this email I realise that I could technically solve 40 
times the equation, only changing the direction each time since the 
components are not coupled to each other (they are howerver coupled with 
another equation not mentionned here) but that looks cumbersome at best and 
probably inefficiency.

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/71014d01-54be-4091-a4b2-594499a5aeb7n%40googlegroups.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

2021-07-26 Thread Sylvain Mathonnière

Thank you for the github link. In the end I fixed the issue by checking 
which value of the dof_index vector was not equal to 
numbers::invalid_unsigned_int. I added this line :

*unsigned int component_i = (dof_index[0] != numbers::invalid_unsigned_int) ? *
*
fe_RTE.system_to_component_index(dof_index[0]).first:*
*
fe_RTE.system_to_component_index(dof_index[1]).first;*

Thank you Corbin as well. I will try that next time I am in need of it !

Best,

Sylvain M.
El viernes, 23 de julio de 2021 a la(s) 20:58:25 UTC+2, Wolfgang Bangerth 
escribió:

> On 7/23/21 9:41 AM, Sylvain Mathonnière wrote:
> > 
> > As for the 4294967295 occurence, is it because I am dealing with an 
> interface 
> > which is a boundary as well ? Isn't this case handled by the 
> *boundary_worker* 
> > and not the *face_worker* (using step-12 notation) ?
>
> No, it's because you are using DG elements. I tried to clarify this here:
> https://github.com/dealii/dealii/pull/12595
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/83b84ffd-3a84-43c6-b918-bb9cc35705bfn%40googlegroups.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

2021-07-23 Thread Sylvain Mathonnière

 Oh, I think I know what those values appears and it is my bad ...
It is probably not 22116 but 22 and 116 and it is probably not 89130 but 
rather 89 and 130. I used something like *std::cout << dof_index[0] << 
std::endl* to get those values and sometimes they just get attached to each 
other like if the *std::endl* was skipped (not sure why, maybe because I 
have several threads ?). I did not realise it before reading your message 
but it actually makes sense because I was not getting an error for those 
values and my Assert was not triggered. I should have seen it, sorry.

As for the 4294967295 occurence, is it because I am dealing with an 
interface which is a boundary as well ? Isn't this case handled by the 
*boundary_worker* and not the *face_worker* (using step-12 notation) ?

Best,

Sylvain

El viernes, 23 de julio de 2021 a la(s) 17:17:17 UTC+2, Wolfgang Bangerth 
escribió:

> On 7/23/21 8:59 AM, Sylvain Mathonnière wrote:
> > Thank you for the answer. I have been looking at a way of doing what you 
> > suggested and I found the function*FEInterfaceValues 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cbe067ec5c8ab4601e70308d94dea8304%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637626491881337201%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000=hrnX9t7udM50eZvFVTiCBunh8ShGkow1iyI4AZq7udw%3D=0>::interface_dof_to_dof_indices()*
>  
>
> > which seems to be doing what I want.
>
> Ah yes, I had forgotten about this function!
>
>
> > I provide it with the interface index "i" and it returns a pair of local 
> > indices corresponding to the active and neighbouring cell. If then I 
> query the 
> > first element *dof_index[0]*, I obtain the local index of the current 
> cell, so 
> > I though of using it like this :
> > 
> > *//get the local indices of the current cell and 
> neighbouring cell*
> > **
> > *std::array dof_index = 
> > fe_iv.interface_dof_to_dof_indices(i);*
> > **
> > **
> > **
> > *   //use local index to obtain component of shape function.*
> > **
> > *   unsigned int component_i = 
> > fe_RTE.system_to_component_index(dof_index[0]).first;*
> > 
> > Then with the debugger and printing some values I realised that 
> dof_index[0] 
> > is kind of always below 160 as expected, but on a couple of occasion it 
> goes 
> > to 22116; 89130 or 4294967295 and triggers the error of 
> > system_to_component_index(). Those values looks like global indices 
> whereas 
> > interface_dof_to_dof_indices(i) should returns local_indices. I must be 
> doing 
> > something horribly wrong but I could not find a tutorial that uses 
> > interface_dof_to_dof_indices(). Am I confused with how to use those 
> functions ?
>
> No, I think you're right that that is how the function is supposed to work.
>
> 4294967295 is numbers::invalid_dof_index, which is the case discussed in 
> the 
> documentation where a DoF lives only one one side of the interface. The 
> other 
> values (22116; 89130) shouldn't happen. If that's what you get, then there 
> is 
> a bug somewhere.
>
> Do you think you could construct a small testcase that illustrates what 
> you have?
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/55622014-8692-448a-92ab-818c134c920an%40googlegroups.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

2021-07-23 Thread Sylvain Mathonnière

 Thank you for the answer. I have been looking at a way of doing what you 
suggested and I found the function* FEInterfaceValues 
::interface_dof_to_dof_indices()*
 
which seems to be doing what I want. 
I provide it with the interface index "i" and it returns a pair of local 
indices corresponding to the active and neighbouring cell. If then I query 
the first element  *dof_index[0]*, I obtain the local index of the current 
cell, so I though of using it like this :

*//get the local indices of the current cell and neighbouring 
cell*
*std::array dof_index = 
fe_iv.interface_dof_to_dof_indices(i);*
 
*   //use local index to obtain component of shape function.*
*   unsigned int component_i = 
fe_RTE.system_to_component_index(dof_index[0]).first;*

Then with the debugger and printing some values I realised that 
dof_index[0] is kind of always below 160 as expected, but on a couple of 
occasion it goes to 22116; 89130 or 4294967295 and triggers the error of 
system_to_component_index(). Those values looks like global indices whereas 
interface_dof_to_dof_indices(i) should returns local_indices. I must be 
doing something horribly wrong but I could not find a tutorial that uses 
interface_dof_to_dof_indices(). Am I confused with how to use those 
functions ?

Best,

Sylvain

El jueves, 22 de julio de 2021 a la(s) 21:53:19 UTC+2, Wolfgang Bangerth 
escribió:

>
> > Here *n_dofs = 320* and so *system_to_component_index()* does not like 
> to have 
> > an argument i>=160 and throws an error.
> > How come n_dofs=320 in this case ? I was expecting 2 interfaces * 2 
> nodes * 40 
> > components = 160 dofs, unless I am considering all the interfaces of the 
> cell 
> > at onces in which case I have 2 interfaces * 4 nodes * 40 components = 
> 320.
> > 
> > *My main question is therefore:*
> > How can I access the ith component of my shape function in face_worker 
> like in 
> > the cell_worker if I cannot use system_to_component_index() ?
>
> You are using FEInterfaceValues, which considers the shape functions that 
> live 
> on the interface. This set of shape functions is the union of the shape 
> functions living on the two adjacent cells. Because you have a 
> discontinuous 
> element, this set has size 2 * the number of dofs per cell, which gives 
> you 
> exactly the 2*160=320 you observe.
>
> Since you can't use FiniteElement::system_to_component_index on your index 
> 'i', you need to find a different way. One approach is to compute a 
> mapping 
> from this index 'i' to a pair '(here_or_there, index within cell)' and 
> then 
> call system_to_component_index() on the second half of the pair.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2dcbdfa0-8a7f-4520-8614-558df5cf2287n%40googlegroups.com.

[deal.II] system_to_component_index() usage in face_worker for DG method

2021-07-22 Thread Sylvain Mathonnière

 Hello everyone !

I have (another) a problem in my code. I am trying to mix *step-12* and 
*step-31* for my project.
I want to solve two coupled equations (named *eq1* and *eq2*) that I am 
solving in parallel just like in *step-31* but for one of them I am using 
the Discontinuous Galerkin method like in *step-12* (*eq1*) and for the 
other one I use standard Galerkin (*eq2*). I also have two *dof_handler* 
like in *step-31* one for each equation. It looks like this :

 


*FESystemfe_eq1; //FE element for equation 1 
FESystemfe_eq2; //FE element for equation 2 
DoFHandler  dof_handler_eq1; DoFHandler  
dof_handler_eq2;*

And the constructor of the class :


*  template  DG_FEM::DG_FEM (): fe_eq1 
(FE_DGQ(1),40), fe_eq2 (FE_Q(1), 1), dof_handler_eq1 
(triangulation), dof_handler_eq2 (triangulation) {}*

For the* eq1*, I am using like in *step-12* the *MeshWorker::mesh_loop(...)* 
to fill my matrix :


In the *cell_worker* function, in order to access the ith component of my 
shape function, I use the same technique as in *step-8*:

...










*const unsigned int n_dofs_per_cell = 
scratch_data.fe_values.get_fe().dofs_per_cell;for (unsigned int i = 
0; i < n_dofs_per_cell; ++i){unsigned int component_i = 
fe_eq1.system_to_component_index(i).first;// and then I access 
the "component_i" of my shape function "i" at gauss point "point" using
 copy_data.cell_rhs(i) +=  fe_v.shape_value_component(i, point, 
component_i) ...}*

...

This seems to work fine. Since I have 40 components and 4 nodes, I have 
*n_dofs_per_cell 
= 160*. That is confirmed by the debuggeur.

My problem arises when I am using the same technique in the *face_worker*.

...









*const unsigned int n_dofs= 
fe_iv.n_current_interface_dofs();for (unsigned int i = 0; i < 
n_dofs; ++i){  unsigned int component_i = 
fe_eq1.system_to_component_index(i).first; 
copy_data_face.cell_matrix(i, j) += * fe_iv.shape_value((beta_dot_n > 0), 
j, qpoint, component_i)}*

...

Here *n_dofs = 320* and so *system_to_component_index()* does not like to 
have an argument i>=160 and throws an error. 
How come n_dofs=320 in this case ? I was expecting 2 interfaces * 2 nodes * 
40 components = 160 dofs, unless I am considering all the interfaces of the 
cell at onces in which case I have 2 interfaces * 4 nodes * 40 components = 
320.

*My main question is therefore:*
How can I access the ith component of my shape function in face_worker like 
in the cell_worker if I cannot use system_to_component_index() ?

Best regards,

Sylvain Mathonnière

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/80dd4f9c-e167-4644-80d0-57a7ee6ae73cn%40googlegroups.com.

Re: [deal.II] How to add 2nd iterator to MeshWorker::mesh_loop ?

2021-07-16 Thread Sylvain Mathonnière

Adding the & did the trick indeed. Thank you. Thank you as well for the FAQ 
entry on this. I completely overlooked it...

Best,

Sylvain

El viernes, 16 de julio de 2021 a la(s) 07:22:25 UTC+2, Jean-Paul Pelteret 
escribió:

> To add to what Wolfgang has already said, there’s this entry in our FAQ on 
> this topic:
>
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#can-i-convert-triangulation-cell-iterators-to-dofhandler-cell-iterators
>
>
> On 15. Jul 2021, at 17:31, Wolfgang Bangerth  
> wrote:
>
> On 7/15/21 2:45 AM, Sylvain Mathonnière wrote:
>
> *typenameDoFHandler::active_cell_iteratorcell_2 
> (cell_1->get_triangulation(), 
> cell_1->level(), cell_1->index(), dof_handler_2);*
>
>
> The compiler error tells you that you need pointers, not references. So 
> this should work:
>
>  typename DoFHandler::active_cell_iterator
> cell_2 (_1->get_triangulation(), cell_1->level(), cell_1->index(),
> _handler_2);
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth  email: bang...@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+un...@googlegroups.com.
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/53c90809-04b1-c2f0-3d35-ec11726fa46f%40colostate.edu
> .
>
>
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d664107d-cc5d-4f26-a543-86c16a6508fdn%40googlegroups.com.

Re: [deal.II] How to add 2nd iterator to MeshWorker::mesh_loop ?

2021-07-15 Thread Sylvain Mathonnière

 I tried as you suggested but it seems there is a little error still. I 
found the constructor in TriaActiveIterator (5/8 constructor) though. 

I wrote :

*typename DoFHandler::active_cell_iterator cell_2 
(cell_1->get_triangulation(), cell_1->level(), cell_1->index(), 
dof_handler_2);*

and the error says (translated from french):
no matching function for the call to *« 
dealii::TriaActiveIterator, 
false> >::TriaActiveIterator(const dealii::Triangulation<2, 2>&, int, int, 
const dealii::DoFHandler<2, 2>&) » 
*
so this apparently match what I wrote.

A bit further down in the error when enumerating the 8 existing constructor 
I got 
*"no known conversion to convert argument 1 from 
« const dealii::Triangulation<2, 2> » to 
« const dealii::Triangulation<2, 2>* »*

So the constructor is expected a pointer (as the doc of TriaActiveIterator 
says) but is getting a value. However I checked the doc of 
*get_triangulation()* and it  return a reference so I am confused now.

Also for the creation of the second iterator, I have to indicate at some 
point the dof_handler I want to link to it, so I guess the last argument 
should be* dof_handler_2* rather than *cell_1-> get_dof_handler()*, 
otherwise there is no difference between iterator cell_1 and cell_2 as far 
as I can tell ? 

Best,

Sylvain

El martes, 13 de julio de 2021 a la(s) 17:26:49 UTC+2, Wolfgang Bangerth 
escribió:

> On 7/13/21 4:02 AM, Sylvain Mathonnière wrote:
> > This sounds like what I need to do indeed. However, I cannot find a 
> function 
> > with such arguments in the help of *DoFHandler*. The closest I could 
> find is 
> > in *DoFAccessor*
> > 
> > *DoFAccessor (
> > const Triangulation< dim, spacedim > *  tria,
> > const int  level,
> > const int  index,
> > const DoFHandler< dim, spacedim > *  dof_handler  )*
> > 
> > Is this the one you were refering to ? but then I do not know how to 
> apply it 
> > to obtain an iterator "cell_2" on the 2nd cell.
>
> The call was this:
> DoFHandler::active_cell_iterator cell_2 (cell_1->get_triangulation(),
> cell_1->level(),
> cell_1->index(),
> cell_1->get_dof_handler());
>
> Note that it's the *iterator* type I'm using here. The class that's behind 
> this is TriaActiveIterator, which indeed has the constructor you need here.
>
> Best
> W>
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a0891ad8-260a-446e-b5a9-087224f6f1e4n%40googlegroups.com.

Re: [deal.II] How to add 2nd iterator to MeshWorker::mesh_loop ?

2021-07-13 Thread Sylvain Mathonnière

This sounds like what I need to do indeed. However, I cannot find a 
function with such arguments in the help of *DoFHandler*. The closest I 
could find is in *DoFAccessor*





*DoFAccessor (const Triangulation< dim, spacedim > *  tria,const int  
level,const int  index,const DoFHandler< dim, spacedim > *  dof_handler  )*

Is this the one you were refering to ? but then I do not know how to apply 
it to obtain an iterator "cell_2" on the 2nd cell.

Something like this does not get me an iterator "cell_2" though:




*DoFAccessor(cell_1->get_triangulation(),cell_1->level(),
 cell_1->index(), dof_handler_2)*

and In the documentation of the class *DoFAccessor*, I could not figure out 
another function to do that.

Maybe I am not using the good class though ?

Best,

Sylvain

El lunes, 12 de julio de 2021 a la(s) 17:00:01 UTC+2, Wolfgang Bangerth 
escribió:

> On 7/12/21 4:18 AM, Sylvain Mathonnière wrote:
> > 
> > _*Question:*_
> > How can I modify the cells worker loop (on dof1) initialisation to have 
> a 
> > second iterator linked with a second dof_handler looping at the same 
> time ? or 
> > how can I add a concurrentely running for loop to the cell_worker with a 
> > different iterator ?
> > 
> > I checked in step 16 and 47 as well as in the doc for MeshWorkers and I 
> did 
> > not find a solution to this problem. I hope the question is clear 
> enough, do 
> > not hesitate to inquires more details if I am missing something 
> important.
>
> Sylvain,
> it's not going to be possible to hack this into MeshWorker::mesh_loop() in 
> an 
> easy way. But what you can do is this: Inside your assembler, whenever you 
> are 
> given a cell iterator pointing to a cell in dof_handler_1, you need to 
> find 
> the corresponding iterator for dof_handler_2. You can do this in the 
> following 
> way (up to typos and misspelled function names):
> DoFHandler::active_cell_iterator cell_2 (cell_1->get_triangulation(),
> cell_1->level(),
> cell_1->index(),
> cell_1->get_dof_handler());
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/3804f06e-92d2-4875-8056-db21553909den%40googlegroups.com.

Re: [deal.II] Solving nonlinear heat equation with radiation

2021-07-05 Thread Sylvain Mathonnière

Thank you for the clear answer, I understand my mistake now, it is much 
simpler like that.

On a side note (irrelevant now), when I was talking about orthonormal 
basis, I was refering to something like \int w_i(x) w_j(x)= 
kronecker_{i,j}, which should be doable and in which basis T_n^3 would look 
like *T_{n}^3 = sum([T_{n}]_i^3 * w_i).*

I was also considering to mimic tutorial 15 with the Newton method if this 
strategy (linearising T^4) was not performing well, but with your reply now 
I will stick to my method.

By the way, Isn't there a mistake in tutorial 15 when it says F'(u_n, δu_n) 
= -F(u_n) (just after "...use a damping parameter αn to get better global 
convergence behavior: "). Shouldn't we read F'(u_n, δu_n) *  δu_n = -F(u_n) 
otherwise it is not homogeneous ?

Best,

Sylvain

El viernes, 2 de julio de 2021 a la(s) 17:49:20 UTC+2, Wolfgang Bangerth 
escribió:

>
> > So I applied the Crank Nicolson discretisation first before  multiplying 
> by 
> > test functions and integrating. But then I am however confused how to 
> handle a 
> > term like *(T_{n}^3 * T_{n+1} , w_i)_omega*
> > (where *(A,B)_omega* is the usual notation for the integral over the 
> spatial 
> > region omega of the product A*B and *w_i* is my test function).
> > 
> > If I write *T_{n}* and *T_{n+1}* as a finite linear combinaison of w_i 
> such 
> > that (*T_{n} = sum_i([T_{n}]_i * w_i)*) then *T_{n}^3* is monstruous 
> unless I 
>
> You're thinking of T_n^3 as a polynomial of particular high degree, but we 
> don't actually care about that. We just need to evaluate it at quadrature 
> points, and to do that you evaluate
> T_n(x_q) = \sum [T_n]_j \varphi_j(x_q)
> and then take whatever value that gives you to the third power. You never 
> need 
> T_n^3 as a *function*, you just need to be able to evaluate it at 
> quadrature 
> points, and so questions such as whether functions are orthogonal don't 
> actually matter.
>
>
> > 1- Is this even mathematically sound ? I should be able to choose my 
> basis 
> > orthogonal as far as I can think. Can I then do that with deal.II ?
>
> Think about it this way: You want to compute an integral
> \int c(x) w_i(x) w_j(x)
> where c(x) is a coefficient that in your case involves T_n^3. You won't in 
> general be able to find orthogonal basis functions w_k for these kind of 
> weights c(x). You just have to approximate the integral by quadrature.
>
>
> > 2 - Is this the approach I should be following for handling such term or 
> > should I rethink the time discretisation scheme ?
>
> No. Just think of the old solution as a coefficient. You probably want to 
> take 
> a look at step-15, which shares many of the characteristics of what you 
> want 
> to do.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/ed8c943c-5df4-49e0-8d19-d9e1f8a049c8n%40googlegroups.com.

Re: [deal.II] Solving nonlinear heat equation with radiation

2021-06-29 Thread Sylvain Mathonnière

I realised I made some errors in my previous post. I need an *orthonormal 
*basis 
of test function (not just orthogonal) to do what I claimed and even then I 
have *T_{n}^3 = sum([T_{n}]_i^3 * w_i) *and then  *(T_{n}^3*T_{n+1}, 
w_j)_omega = sum_l( T_{n+1}_l (sum_k( T_{n}_k^3 * w_k*w_l, w_j))_omega).* 

El martes, 29 de junio de 2021 a la(s) 13:24:31 UTC+2, Sylvain Mathonnière 
escribió:

> Thank you for your answer ! I was trying to follow the paper where they do 
> spatial discretisation before time discretisation, so that is why.
> Why would one typically prefer to do it the other way (time discretisation 
> before spatial) ? I checked and it is indeed the case in the deal.II 
> tutorials I looked at.
>
> I realised now that my problem is not completely linked with deal.II and 
> so my questions might seem inapproriate on such forum but I saw in the 
> guideline that this forum is " is also open to some broader discussions 
> on the finite element method, numerical methods, C++ and similar". So I 
> will give my questions a shot ^^. (don't hesitate to politely send me back 
> to another place if this is getting too far for the purpose of this forum).
>
> So I applied the Crank Nicolson discretisation first before  multiplying 
> by test functions and integrating. But then I am however confused how to 
> handle a term like *(T_{n}^3 * T_{n+1} , w_i)_omega*
> (where *(A,B)_omega* is the usual notation for the integral over the 
> spatial region omega of the product A*B and *w_i* is my test function).
>
> If I write *T_{n}* and *T_{n+1}* as a finite linear combinaison of w_i 
> such that (*T_{n} = sum_i([T_{n}]_i * w_i)*) then *T_{n}^3* is monstruous 
> unless I used an orthogonal basis of test function *(w_i,w_j) 
> =delta_{i,j}* and then I end up with *T_{n}^3 = sum([T_{n}]_i^3 * w_i)*. 
> If I follow my orthogonal basis for test function (is it even legal ?) I 
> end up with:
>
> for 1<= j <= total_dof
> *(T_{n}^3*T_{n+1}, w_j)_omega = sum_l( T_{n+1}_l (sum_k( T_{n}_k * 
> w_k^3*w_l, w_j))_omega)*, which is not as bad. So I would have a system 
> like A*T_{n+1} where A is a matrix with element (l,j) : *sum_k( T_{n}_k * 
> w_k^3*w_l, w_j))_omega.*
>
> If I follow this method, I end up with a linear system. *However I then 
> have 2 questions :*
>
> 1- Is this even mathematically sound ? I should be able to choose my basis 
> orthogonal as far as I can think. Can I then do that with deal.II ?
> 2 - Is this the approach I should be following for handling such term or 
> should I rethink the time discretisation scheme ?
>
> Finally, I want to thank you for your time as I notice you personally 
> answer many questions on this forum and many others in other forums (I 
> found you in two more forums !)
>
> Best,
>
> Sylvain Mathonnière
> El martes, 29 de junio de 2021 a la(s) 03:58:58 UTC+2, Wolfgang Bangerth 
> escribió:
>
>>
>> Sylvain, 
>> I have not tried to follow your equations in detail, but let me point out 
>> how 
>> this is generally done. The "right" approach is to start with the PDE, 
>> then to 
>> use a time discretization *of the PDE* that leads to a linear occurrence 
>> of 
>> T^{n+1} on the left hand side (possibly multiplied by factors that 
>> involve 
>> T^n), and only then think about the spatial discretization. For the 
>> latter, as 
>> usual you multiply by a test function and integrate by parts, then 
>> approximate 
>> the integrals by quadrature. 
>>
>> From your description, it seems to me that you are trying to force a 
>> special 
>> structure of your linear system, and what you arrive at may or may not be 
>> correct, but you don't know why because you are skipping steps of the 
>> outline 
>> above. For example, if you have a term of the form 
>>
>> T_n^3 T_{n+1} 
>>
>> then in general this will *not* lead to multiplying the vector of 
>> unknowns 
>> \vec T_{n+1} by a diagonal matrix with entries equal to [T_n]_i^3 *unless 
>> you 
>> use a particular quadrature formula*. Rather, if you use a Gauss 
>> quadrature 
>> formula, you will end up multiplying [T_{n+1}]_i by factors that will 
>> also 
>> include [T_n]_j for degrees of freedom j that are neighbors of i. 
>>
>> Best 
>> W. 
>>
>>
>> On 6/28/21 8:11 AM, Sylvain Mathonnière wrote: 
>> > 
>> > Dear all, 
>> > 
>> > *Some background :* 
>> > ** 
>> > I have been working for some time now on a project which is to solve 
>> the 
>> > *radiative transfer equation* using the discrete ordinate method 
>> coupled with 
>>

Re: [deal.II] Solving nonlinear heat equation with radiation

2021-06-29 Thread Sylvain Mathonnière

 Thank you for your answer ! I was trying to follow the paper where they do 
spatial discretisation before time discretisation, so that is why.
Why would one typically prefer to do it the other way (time discretisation 
before spatial) ? I checked and it is indeed the case in the deal.II 
tutorials I looked at.

I realised now that my problem is not completely linked with deal.II and so 
my questions might seem inapproriate on such forum but I saw in the 
guideline that this forum is " is also open to some broader discussions on 
the finite element method, numerical methods, C++ and similar". So I will 
give my questions a shot ^^. (don't hesitate to politely send me back to 
another place if this is getting too far for the purpose of this forum).

So I applied the Crank Nicolson discretisation first before  multiplying by 
test functions and integrating. But then I am however confused how to 
handle a term like *(T_{n}^3 * T_{n+1} , w_i)_omega*
(where *(A,B)_omega* is the usual notation for the integral over the 
spatial region omega of the product A*B and *w_i* is my test function).

If I write *T_{n}* and *T_{n+1}* as a finite linear combinaison of w_i such 
that (*T_{n} = sum_i([T_{n}]_i * w_i)*) then *T_{n}^3* is monstruous unless 
I used an orthogonal basis of test function *(w_i,w_j) =delta_{i,j}* and 
then I end up with *T_{n}^3 = sum([T_{n}]_i^3 * w_i)*. If I follow my 
orthogonal basis for test function (is it even legal ?) I end up with:

for 1<= j <= total_dof
*(T_{n}^3*T_{n+1}, w_j)_omega = sum_l( T_{n+1}_l (sum_k( T_{n}_k * 
w_k^3*w_l, w_j))_omega)*, which is not as bad. So I would have a system 
like A*T_{n+1} where A is a matrix with element (l,j) : *sum_k( T_{n}_k * 
w_k^3*w_l, w_j))_omega.*

If I follow this method, I end up with a linear system. *However I then 
have 2 questions :*

1- Is this even mathematically sound ? I should be able to choose my basis 
orthogonal as far as I can think. Can I then do that with deal.II ?
2 - Is this the approach I should be following for handling such term or 
should I rethink the time discretisation scheme ?

Finally, I want to thank you for your time as I notice you personally 
answer many questions on this forum and many others in other forums (I 
found you in two more forums !)

Best,

Sylvain Mathonnière
El martes, 29 de junio de 2021 a la(s) 03:58:58 UTC+2, Wolfgang Bangerth 
escribió:

>
> Sylvain,
> I have not tried to follow your equations in detail, but let me point out 
> how 
> this is generally done. The "right" approach is to start with the PDE, 
> then to 
> use a time discretization *of the PDE* that leads to a linear occurrence 
> of 
> T^{n+1} on the left hand side (possibly multiplied by factors that involve 
> T^n), and only then think about the spatial discretization. For the 
> latter, as 
> usual you multiply by a test function and integrate by parts, then 
> approximate 
> the integrals by quadrature.
>
> From your description, it seems to me that you are trying to force a 
> special 
> structure of your linear system, and what you arrive at may or may not be 
> correct, but you don't know why because you are skipping steps of the 
> outline 
> above. For example, if you have a term of the form
>
> T_n^3 T_{n+1}
>
> then in general this will *not* lead to multiplying the vector of unknowns 
> \vec T_{n+1} by a diagonal matrix with entries equal to [T_n]_i^3 *unless 
> you 
> use a particular quadrature formula*. Rather, if you use a Gauss 
> quadrature 
> formula, you will end up multiplying [T_{n+1}]_i by factors that will also 
> include [T_n]_j for degrees of freedom j that are neighbors of i.
>
> Best
> W.
>
>
> On 6/28/21 8:11 AM, Sylvain Mathonnière wrote:
> > 
> > Dear all,
> > 
> > *Some background :*
> > **
> > I have been working for some time now on a project which is to solve the 
> > *radiative transfer equation* using the discrete ordinate method coupled 
> with 
> > the *nonlinear heat equation* which include radiative terms. I am 
> following 
> > this approach here
> > 
> > https://hal.archives-ouvertes.fr/hal-01273062/document
> > 
> > At each time step, I am solving first the radiative transfer equation 
> and then 
> > solving the heat equation and then moving to the next time step.
> > 
> > I am mainly following tutorial 31 since it seems close to my problem 
> > (replacing Navier-Stokes with Radiative transfer equation and adding the 
> > nonlinear term in the heat equation), I used it to see how to couple 
> equations 
> > in the solver and how to obtain result for previous time step.
> > 
> > Among the particularities of my program is the use of discontinuous 
> galerkin 
> > method for the radiative transfer equation (RT

[deal.II] Re: Solving nonlinear heat equation with radiation

2021-06-28 Thread Sylvain Mathonnière

It seems my image of the 2D version did not get through. It was something 
like  (matrix * vector = vector):

| T_{n}_10  |  |T_{n+1}_1|   =   |T_{n}_1 * T_{n+1}_1|
| 0   T_{n}_2   |  |T_{n+1}_2||T_{n}_2 * T_{n+1}_2 |

I hope it goes through nicely...


-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6329fc29-2337-4c34-a34b-e2091d23a449n%40googlegroups.com.

[deal.II] Solving nonlinear heat equation with radiation

2021-06-28 Thread Sylvain Mathonnière


Dear all,

*Some background :*

I have been working for some time now on a project which is to solve the 
*radiative 
transfer equation* using the discrete ordinate method coupled with the 
*nonlinear 
heat equation* which include radiative terms. I am following this approach 
here 

https://hal.archives-ouvertes.fr/hal-01273062/document

At each time step, I am solving first the radiative transfer equation and 
then solving the heat equation and then moving to the next time step. 

I am mainly following tutorial 31 since it seems close to my problem 
(replacing Navier-Stokes with Radiative transfer equation and adding the 
nonlinear term in the heat equation), I used it to see how to couple 
equations in the solver and how to obtain result for previous time step.

Among the particularities of my program is the use of discontinuous 
galerkin method for the radiative transfer equation (RTE) and more standard 
galerkin method for the heat equation (HE), my constructor looks like this :

template  DG_FEM::DG_FEM ()
:fe_HE (FE_Q(1), 1), dof_handler_HE (triangulation),
fe_RTE (FE_DGQ(1), N_dir), dof_handler_RTE (triangulation)
{}
With fe_HE and fe_RTE as FEsystem. Here N_dir represents the direction 
discretisation used in the discrete ordinate method and for my case I have 
N_dir = 40. This is the best way I could come up to for this problem.

I use two dof_handler since I uses FE_Q and FE_DGQ like in tutorial 31 and 
so it is not too confusing for me.

I hope this part is relatively clear, if not, or if it is clear but 
something looks horrifiyingly wrong with the way I am doing, please do not 
hesitate to inquire more from me.


*Getting to my question :*

When time comes to solve the heat equation I use the Crank-Nicolson scheme 
(like in the paper, starting equation 44), but at the end, I have trouble 
to obtain a linear system like AT = B (A=total matrix, T=temperature, B= 
RHS). The paper is really vague on this part and I believe the notation is 
misleading from equation 56.

The way it is done is by linearising the nonlinear term T^4 using Taylor 
Young approximation : 
T_{n+1}^4 = T_{n}^4 + 4 T_{n}^3 * (T_{n+1} - T_{n})
Keep in mind T_{n+1} are vectors and the "*" represents term to term 
multiplication.

Following the article I end up with the heat equation looking like :

[ M/k + A/2] T_{n+1} = [ M/k - A/2] T_{n} - theta * M * (-T_{n}^4 + 2 **T_{n}^3 
* T_{n+1}*) + extra_RHS

Where M is the mass matrix, A is the stiffness matrix, k the time step and 
everything else is constant not relevant to my question.

The annoying part of this equation is that I cannot factor all the term 
T_{n+1} easily on the left side since "T_{n}^3 * T_{n+1}" is actually a 
vector (T_{n}_i^3*T_{n+1}_i) for i €[1;total_dof]. The only way I found 
(and I hope it is correct) is to rewrite :

*(T_{n}^3 * T_{n+1})_{1<=i<=total_dof} = Diag(T_{n}^3) * T_{n+1}.*
A 2D version of what I would like to do is shown below:


I am basically creating a diagonal matrix with (T_{n}^3)_i on the diagonal; 
The product of this diagonal matrix with the vector T_{n+1} should give me 
back the vector I want.
If I can do that then I can factor the vector T_{n+1} and rewrite my 
equation as :

[ M/k + A/2 + 2 * Diag(T_{n}^3) ] T_{n+1} = Function_of( T_{n}) 

and I am happy.

*My question is :* How can I do this Diagonal matrix creation in a setting 
similar to tutorial 31 where I am looping through the cells and filling a 
local_matrix everytime (so I do not have access to the full T_{n} vector 
but just the local cell T_{n} vector) ?

Finally, generally speaking, does this approach sounds "healthy" to you ?  
or should I use some Newton's approach like in tutorial 15 ? or is there 
something obvious I am missing, like my math is wrong or something else ?

Any help would be greatly appreciated.

Best regards,

Sylvain Mathonnière


-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/bdcb220c-251a-486d-a047-b7030597383dn%40googlegroups.com.

Re: [deal.II] integrating over boundaries

2021-04-26 Thread Sylvain Mathonnière

Thank you very much for your help ! Indeed it is simpler outside of the 
DataPostprocessor class. 
I compared the integration using quadrature points and using constant value 
along cell face (as I initially intended) and got similar results which 
comfort me that I have the correct implementation.
The key enabler was really to get the gradient of the solution vector using 
"fe_face_values.get_function_gradients(solution_vector, gradients)" for me. 
Thanks again for the help.

Best regards,

Sylvain Mathonnière

El domingo, 25 de abril de 2021 a la(s) 20:50:05 UTC+2, Jean-Paul Pelteret 
escribió:

> Hi Sylvain,
>
> Thanks for detailing your attempts so thoroughly. With this
>
>   // I recover the value of conductivity
>   MaterialData material_data;
>
>   const typename DoFHandler::cell_iterator current_cell = 
> input_data.template get_cell>();
>
>   FullMatrix  conductivity = 
> material_data.get_th_conductivity(current_cell->material_id());
>
>   double heat_flux_cell_norm = 0.;
>
>   for (auto face_index : GeometryInfo::face_indices())
>   {
>   if (current_cell->face(face_index)->boundary_id() == 2)
>   {
> double long_cell = 0;
>
> long_cell = (current_cell->face(face_index)->vertex(0) - 
> current_cell->face(face_index)->vertex(1)).norm();
>
> for (unsigned int d=0; d {
>
>   heat_flux_cell_norm += 
> -conductivity[d][d]*input_data.solution_gradients[face_index][d] * 
>
>   
> conductivity[d][d]*input_data.solution_gradients[face_index][d];
> }
> 
> std::cout << "total heat flux = " << heat_flux_cell_norm << std::endl;
> heat_flux_cell_norm *= long_cell;
>   }
>   }
>
>
> you’ve got most of the structure for your calculation completely correct, 
> but embedding it within a DataPostprocessor is not what you want to do. 
> These classes are there to help visualise quantities via DataOut (i.e. do 
> some per-cell calculations), rather than computing some boundary integral 
> like you want. Instead, you should encapsulate the above within your own 
> function and the standard DoFHandler active cell loop. You would then 
> replace “input_data” with some FEFaceValues object (reinitialised on the 
> given cell and face_index) and then call "
> fe_face_values.get_function_gradients 
> <https://dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#ad1f4e0deb5d982e8172d82141c634a67>(solution_vector,
>  
> gradients);” to store the solution gradients at each face quadrature point 
> and from that compute the heat flux at each face quadrature point. But I 
> think that you’re also just wanting to consider the component of the heat 
> flux that is normal to the boundary, right? 
>
> I think that your calculation should have this sort of structure:
>
> double total_heat_flux = 0.0;
> FEFaceValues fe_face_values(fe, face_quadrature, update_gradients | 
> update_normal_vectors | update_JxW_values);
> std::vector> temperature_gradient(face_quadrature.size());
>
>   for (const auto cell : dof_handler.active_cell_iterators())
> for (const auto face_index : GeometryInfo::face_indices())
>   If (… at boundary of interest …)
>   {
>fe_face_values.reinit(cell, face_index);
>fe_face_values.reinit.get_solution_gradents(solution_vector, 
> temperature_gradient);
> for (const auto q_point : 
> fe_face_values.quadrature_point_indices())
> {
>   const Tensor<1,dim> heat_flux_q = 
> -conductivity*temperature_gradient[q_point]; // Q = -K.Grad(T)
>   total_heat_flux += heat_flux_q 
> * fe_face_values.normal_vector(q_point) * fe_face_values.JxW(q_point); // q 
> = Q.N. ; q_total = \int q dA
>   }
> }
>
> If I’ve interpreted things correctly then that (or something close to it) 
> should achieve what you’re looking to do.
>
> I hope that this help!
>
> Best,
> Jean-Paul
>
> On 23. Apr 2021, at 12:03, Sylvain Mathonnière  
> wrote:
>
> By the way, I am aware than I forgot to take the square root of 
> heat_flux_cell_norm. I am also aware that using 
> *input_data.solution_gradients[face_index][d]*  is weird in this context 
> since the first indices should refer to a vertex indices of the cell and 
> not the face index and it only works because quadrangle have 4 faces and 4 
> vertices but so far this was not my main issue but it is definitively on my 
> list of fixing.
>
> El viernes, 23 de abril de 2021 a la(s) 10:07:23 UTC+2, Sylvain 
> Mathonnière escribió:
>
>> Hello ever

[deal.II] integrating over boundaries

2021-04-23 Thread Sylvain Mathonnière

Hello everyone !

I have a very basic question but somehow I did not manage to find the 
answer in tutorials, neither looking at the documentation or in the this 
user group. Or more probably I did not understand it when I saw it since I 
am fairly new to deal.II and to C++.

*My main goal :* I am trying to integrate over a boundary (a line since I 
am in 2D) the heat flux to get the total heat flux through my boundary.

*What I did :*
1 - I looked at this documentation 
https://www.dealii.org/current/doxygen/deal.II/classDataPostprocessorVector.html
 
to calculate first the heat flux on my boundary which worked fine.

2 - I tried to integrate the norm of the heat flux over my boundary in the 
function *evaluate_scalar_field()* function but I could not. My best effort 
is this :

template 
class HeatFluxPostprocessor : public DataPostprocessorVector
{
public:
  HeatFluxPostprocessor ()
:
DataPostprocessorVector ("_k_grad_T",
  update_gradients | update_quadrature_points)
  {}
  virtual void evaluate_scalar_field(
const DataPostprocessorInputs::Scalar _data,
std::vector >   _quantities) const
{
  // I recover the value of conductivity
  MaterialData material_data;
  const typename DoFHandler::cell_iterator current_cell = 
input_data.template get_cell>();
  FullMatrix  conductivity = material_data.get_th_conductivity(
current_cell->material_id());

  double heat_flux_cell_norm = 0.;

  for (auto face_index : GeometryInfo::face_indices())
  {
  if (current_cell->face(face_index)->boundary_id() == 2)
  {
double long_cell = 0;
long_cell = (current_cell->face(face_index)->vertex(0) - 
current_cell->face(face_index)->vertex(1)).norm();

for (unsigned int d=0; dhttp://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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6058222b-5470-4bbd-ad85-7f43327b37adn%40googlegroups.com.

Re: [deal.II] Loss of symmetry by changing parameter

Re: [deal.II] Loss of symmetry by changing parameter

[deal.II] Re: Loss of symmetry by changing parameter

Re: [deal.II] Re: function face_component_to_system_index() for Discontinuous Galerkin ?

[deal.II] Re: function face_component_to_system_index() for Discontinuous Galerkin ?

[deal.II] function face_component_to_system_index() for Discontinuous Galerkin ?

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[deal.II] system_to_component_index() usage in face_worker for DG method

Re: [deal.II] How to add 2nd iterator to MeshWorker::mesh_loop ?

Re: [deal.II] How to add 2nd iterator to MeshWorker::mesh_loop ?

Re: [deal.II] How to add 2nd iterator to MeshWorker::mesh_loop ?

Re: [deal.II] Solving nonlinear heat equation with radiation

Re: [deal.II] Solving nonlinear heat equation with radiation

Re: [deal.II] Solving nonlinear heat equation with radiation

[deal.II] Re: Solving nonlinear heat equation with radiation

[deal.II] Solving nonlinear heat equation with radiation

Re: [deal.II] integrating over boundaries

[deal.II] integrating over boundaries

20 matches

Site Navigation

Mail list logo

Footer information