Different shape representations with manifolds on the same triangulation

[Juan Carlos Araujo Cabarcas]

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

[Jean-Paul Pelteret]

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 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() function could be helpful for this purpose. This, though, seems less elegant than the first approach.

Best,

Jean-Paul

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/92fc053d-ca4e-4830-bca0-2f211a9ba384%40googlegroups.com.

[Juan Carlos Araujo Cabarcas]

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,

To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[Jean-Paul Pelteret]

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

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

To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/e1b5d8e7-97cb-4bf8-a562-1d7aa9429536%40googlegroups.com.

[Juan Carlos Araujo Cabarcas]

Dear Jean-Paul, thanks again for your support and kind suggestions.

I have worked with MappingQEulerian some time before, and as I remember, I dropped it because it requires the use of a vector field defined in the whole domain in order to curve geometries.

After this, I started using ChartManifolds and Transfinite Interpolation, because it is convenient for me to define the geometry only on edges of a primitive, and then bend them.

However, I will look again into MappingQEulerian, maybe there is something of this class that I can use.

Regarding your observations on efficiency, I must admit that I generalized my setting way too much. My functions f(x), g(x) are just constants f0,g0 outside, and f1,g1 inside of the inner shape.

Additionally, the trace/rope of the curved shape coincides with the edges of my elements. Then, evaluating f(x), g(x) is just knowing if the current cell is inside or outside the inner shape, which does not require a complicated transformation.

In my setting, I use FE to approximate the solution u(q) of my PDE, and then I evaluate a given functional J[q,u(q)], for a vector

q=(alpha_1,alpha_2,...,alpha_m,...,alpha_N)^T with geometrical parameters alpha_m, and unitary vectors e_m, m=1,2,....N.

Then my aim is to test with Finite Differences my derivation of the Shape Derivative given in Ch. 9, Sec. 3, and examples in Sec. 4. of the reference:

    Shapes and Geometries. Metrics, Analysis, Differential Calculus, and Optimization. M. C. Delfour, J.-P. Zolesio, 2nd edition, SIAM, 2011

Approximating the Shape derivative implies that , I need to perform the assembly of the mass and stiffness matrices for J(q+eps*e_m), J(q), with 0<eps a small number.

I hope this clarifies my aim with this thread.

Again, it would be nice if I could define N mappings, and at every cell, I could just apply the desired mapping and extract a FEValues object.

I will continue to dig into this, but please feel free to let me know if there is a simpler way.

Maybe, I am complicating things too much.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/e1b5d8e7-97cb-4bf8-a562-1d7aa9429536%40googlegroups.com.