Dear Juan Carlos, > If I interpret your mail correctly, your suggestion is using FunctionManifold > in > order to create the parametrization of the geometry. If this is the case, I > have > already done this step through my own class PolarShapeManifold (no code > included). > If you meant differently, can you please clarify?
I think that you’ve done what I was thinking was necessary, but now that I’ve thought about it some more its the MappingQEulerian <https://dealii.org/current/doxygen/deal.II/classMappingQEulerian_1_1MappingQEulerianGeneric.html> class that provides some mechanism to push forward from a reference geometry to one that the computations are done on. If that is a possibility then you could, conceptually, use (or abuse) the hp framework to have only one triangulation and multiple mappings, and just switch between them. I guess that would work to minimise the number of triangulations that you would need to store. Note again that I’m really just throwing out ideas here, because I’m not too familiar with the manifold classes beyond the generic ones that deal.II provides. I’ve not got a significant amount of familiarity with the mapping classes outside of their standard applications. > I perform N independent assembly routines (loop on the cells of each > triangulation). Ok, well what I think is important to recognise is that certain components of the assembly operation depend on how the real cell gets mapped to the reference element, and some not. This is all sketched out here: https://dealii.org/current/doxygen/deal.II/group__FE__vs__Mapping__vs__FEValues.html <https://dealii.org/current/doxygen/deal.II/group__FE__vs__Mapping__vs__FEValues.html> So, for example, the shape functions will remain unchanged for the I’th DoF evaluated at the q’th quadrature point of each cell, irrespective of real cell geometry. However, the gradients, mapping Jacobian (and, consequently, weighted Jacobian for numerical integration) will not. So, in theory, I guess it might be able to (for example) multiplicatively decompose your assembled matrix into some parts that are geometry independent and others that are, I’m not sure if there’s much to be gained by doing so. For your example (stated in your original post) both contributions to your stiffness matrix appear to be either position or cell-geometry dependent. Additionally, working with such data structures seems so much more inefficient than just assembling the problem in the first place. If you have some geometry-invariant data then, of course, there’s nothing stopping you from pre-computing it and just using it as would be required during the assembly process. In summary, I can’t think of anything useful to help reduce the amount of computation effort during the assembly process for the example you previously described. I guess the logical question for me to ask is, why do you want to do this in the first place? Is the assembly definitely a bottle-neck for you? Best, Jean-Paul > On 16 Jan 2020, at 13:19, Juan Carlos Araujo Cabarcas <ju4...@gmail.com> > wrote: > > Dear Jean-Paul, thank for your interest in my problem and your quick reply! > > If I interpret your mail correctly, your suggestion is using FunctionManifold > in > order to create the parametrization of the geometry. If this is the case, I > have > already done this step through my own class PolarShapeManifold (no code > included). > If you meant differently, can you please clarify? > > > What I currently do, is creating N triangulations, and to each one I apply a > different PolarShapeManifold that gives a shape according to the parameter > alpha_m, m=1,2,...N. > > Consecutively, I perform N independent assembly routines (loop on the cells > of each triangulation). > > As you may notice, I would like to avoid performing redundant operations and > using > extra instantiations (triangulation). > > Then, since all triangulations have the same structure (nodes, edges, > coloring) > it would be convenient to loop only through the cells of a "master" > triangulation and > apply different mappings (alpha_m) to the cell iterator as we loop. > > Maybe this is not currently possible? > > Thanks in advance, > > Juan Carlos Araújo, > > > On Wednesday, 15 January 2020 22:19:46 UTC+1, Jean-Paul Pelteret wrote: > Dear Juan Carlos, > > This is not my area of expertise, so I’m sticking my neck out to offer some > suggestions (hoping that I’ve understood the problem in the first place). > > Might it be possible to do achieve this though use of a customised Manifold > class? There is the FunctionManifold > <https://www.dealii.org/current/doxygen/deal.II/classFunctionManifold.html> > class that seems (to me, at least) to offer the functionality that you’re > requiring. So you wouldn’t modify your assembly code or geometry at all, but > only define this parametric FunctionManifold which then, though the mapping, > morphs the interpretation of the geometry as is seen by FEValues you require. > Naturally, your mapping functions need to be well defined everywhere on the > domain. > > Alternatively, could you simply deform the grid itself as part of a > pre-processing step? The GridTools::laplace_transform() > <https://www.dealii.org/current/doxygen/deal.II/namespaceGridTools.html#a7ed2aaa1aea3ac22b1e1807ce6d0b5f3> > function could be helpful for this purpose. This, though, seems less elegant > than the first approach. > > Best, > Jean-Paul > >> On 15 Jan 2020, at 14:46, Juan Carlos Araujo Cabarcas <ju4...@gmail.com >> <javascript:>> wrote: >> >> Dear all, >> >> I would like your guidance on how to perform the assembly of different >> shape representations on the same triangulation and on the same loop through >> cells. >> >> Let me try to explain a bit more. >> >> I have designed a grid containing two concentric squares (or circles). >> Additionally, I have a parametric representation T(alpha) that maps points >> (x,y) on the edges >> of the inner square (or circle) to points (x',y') on the edges of a desired >> shape >> contained inside the outer square (or circle). This is (x',y')=T(alpha;x,y), >> where >> the representation T depends on a parameter alpha. >> >> Now assume we have N different shapes that are labelled with the index m. >> We then have N parameters alpha_m, with m=1,2,...N. >> Then in the loop through the triangulation cells in the assembly process, I >> would like >> to be able for each cell to loop through different alpha_m, in order to >> generate the >> local FE matrices and load vectors corresponding to each of the >> modified triangulations corresponding to each of the N shapes. >> >> >> My guess is that somehow I should modify the following line: >> >> const FEValues<dim> &fe_v = hp_fe_v.get_present_fe_values (); >> >> >> and pass instead a FEValues that has been generated with the desired alpha_m. >> >> Below I provide a more complete sketch on how my code looks like. >> >> I am quite hesitant on where to start, and I would appreciate your insights >> on >> how to achieve having an assemly routine where I can pass different alpha_m >> per each >> cell in the loop through cells. >> I am grateful for any advise or hint on how to achieve this. >> >> Thanks in advance, >> Juan Carlos Araújo, PhD >> >> >> >> The way I work with one shape is the following: >> >> // Deaclare my environment >> PolarShapeManifold manifold; >> PolarManifold<dim> polar; >> TransfiniteInterpolationManifold<dim> inner_manifold; >> >> Triangulation<dim> triangulation; >> >> hp::DoFHandler<dim> dof_handler; >> hp::FECollection<dim> fe_collection; >> hp::MappingCollection<dim> mapping_collection; >> >> PETScWrappers::SparseMatrix system_matrix; >> PETScWrappers::MPI::Vector solution; >> >> // Constructor >> manifold ( alpha), dof_handler (triangulation),... >> >> // setup_system >> typename hp::DoFHandler<dim>::active_cell_iterator >> cell = dof_handler.begin_active(), >> endc = dof_handler.end(); >> for (; cell!=endc; ++cell) { >> cell->set_active_fe_index ( ... ); >> } >> dof_handler.distribute_dofs (fe_collection); >> ... >> >> // create_coarse_grid >> concentric_disks_inner_shape ( triangulation ); >> tria.set_manifold ( manifold_label, manifold ); >> tria.set_manifold ( layer_label, polar ); >> >> unsigned int not_curved = 1; >> inner_manifold.initialize(triangulation); >> tria.set_manifold (not_curved, inner_manifold); >> ... >> >> // Assembly >> >> hp::FEValues<dim> hp_fe_v ( mapping_collection, >> >> fe_collection, >> quadrature_collection, >> update_values | update_gradients | >> update_quadrature_points | >> update_JxW_values); >> >> FullMatrix<PetscScalar> cell_system; >> Vector<PetscScalar> cell_rhs; >> >> std::vector<types::global_dof_index> local_dof_indices; >> >> typename hp::DoFHandler<dim>::active_cell_iterator >> cell = dof_handler.begin_active(), >> endc = dof_handler.end(); >> for (; cell!=endc; ++cell) { >> const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; >> >> cell_system.reinit (dofs_per_cell, dofs_per_cell); >> cell_rhs.reinit (dofs_per_cell); >> ... >> >> hp_fe_v.reinit (cell); >> const FEValues<dim> &fe_v = hp_fe_v.get_present_fe_values (); >> >> for (unsigned int q_point=0; >> q_point<fe_v.n_quadrature_points; ++q_point) { >> Point<dim> x = fe_v.quadrature_point (q_point); >> >> for (unsigned int i=0; i<dofs_per_cell; ++i){ >> for (unsigned int j=0; j<dofs_per_cell; ++j) { >> >> cell_system (i,j) += ( >> >> fe_v.shape_grad(i,q_point) * >> fe_v.shape_grad(j,q_point) >> -f(x) * // given f >> fe_v.shape_value(i,q_point) * >> fe_v.shape_value(j,q_point) ) * >> fe_v.JxW(q_point); >> } >> >> cell_rhs (i) += g(x) * // given g >> >> fe_v.shape_value(i,q_point) * >> fe_v.JxW(q_point); >> } >> } // q >> >> >> >> >> >> >> >> >> >> -- >> The deal.II project is located at http://www.dealii.org/ >> <http://www.dealii.org/> >> For mailing list/forum options, see >> https://groups.google.com/d/forum/dealii?hl=en >> <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 dea...@googlegroups.com <javascript:>. >> To view this discussion on the web visit >> https://groups.google.com/d/msgid/dealii/92fc053d-ca4e-4830-bca0-2f211a9ba384%40googlegroups.com >> >> <https://groups.google.com/d/msgid/dealii/92fc053d-ca4e-4830-bca0-2f211a9ba384%40googlegroups.com?utm_medium=email&utm_source=footer>. > > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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 > <mailto:dealii+unsubscr...@googlegroups.com>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/e1b5d8e7-97cb-4bf8-a562-1d7aa9429536%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/e1b5d8e7-97cb-4bf8-a562-1d7aa9429536%40googlegroups.com?utm_medium=email&utm_source=footer>. -- 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/798ACFD3-54E4-40E9-9ACF-25D9517E247D%40gmail.com.
