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/71014d01-54be-4091-a4b2-594499a5aeb7n%40googlegroups.com.
