*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<dim> * = *fe_eq1 (FE_DGQ<dim>(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<dim>::active_cell_iterator cell_2 (&cell->get_triangulation(), cell->level(), cell->index(), &dof_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: [email protected] > 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1bc6bd86-3d8e-473f-bbe2-eb8e1873dcf1n%40googlegroups.com.
