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/ > 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 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/ 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/e1b5d8e7-97cb-4bf8-a562-1d7aa9429536%40googlegroups.com.
