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 
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 

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 
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;
    tria.set_manifold (not_curved, inner_manifold);

// Assembly

        hp::FEValues<dim> hp_fe_v ( mapping_collection,
                                update_values    |  update_gradients |
                                update_quadrature_points  | 
    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) += (
shape_grad(i,q_point) *
                                 -f(x) *     // given f
                                 fe_v.shape_value(i,q_point) *
                                 fe_v.shape_value(j,q_point) ) *
                            cell_rhs (i)   +=  g(x) *      // given g
shape_value(i,q_point) *
                    } // q

The deal.II project is located at
For mailing list/forum options, see
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 view this discussion on the web visit

Reply via email to