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<dim>(1)*.
More generally, it should be something like *n_nodes =
(fe_eq1.degree+1)^dim*, so n_nodes = 4 for *FE_DGQ<dim>(1)* but n_nodes =
9 for *FE_DGQ<dim>(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 :
>
> * FESystem<dim> fe_eq1; *
> * DoFHandler<dim> dof_handler_eq1;*
>
> And the constructor of the class :
>
> * template <int dim> DG_FEM<dim>::DG_FEM ()*
> * : fe_eq1 (FE_DGQ<dim>(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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/2334bf8a-a13a-4ec2-875f-333974dd6f59n%40googlegroups.com.
