Re: [deal.II] Different shape representations with manifolds on the same triangulation

[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
> 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
>
> On 16 Jan 2020, at 13:19, Juan Carlos Araujo Cabarcas  > 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 

Re: [deal.II] Different shape representations with manifolds on the same triangulation

[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 a vector field define in the whole domain in 
order to curve the geometry.
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, it is just knowing if the current cell is inside or 
outside the inner shape, which does not require a complicated 
transformation.

Finally, I am testing 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 for a vector 
x=(alpha_1,alpha_2,...,alpha_m,...,alpha_N)^T with geometrical parameters 
alpha_m, and unitary vectors e_m, m=1,2,N, I need to perform the 
assembly for f(x+eps*e_m), f(x).
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.


On Sunday, 19 January 2020 23:50:37 UTC+1, Jean-Paul Pelteret wrote:
>
> 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
>
> On 16 Jan 2020, at 13:19, Juan Carlos Araujo Cabarcas  > 

Re: [deal.II] Different shape representations with manifolds on the same triangulation

[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 a vector field define in the whole domain in 
order to curve the geometry.
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 
there that I could 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, it is just knowing if the current cell is inside or 
outside the inner shape, which does not require a complicated 
transformation.

Finally, I am testing with Finite Differences my derivation of the Shape 
Derivative given in Ch. 9, Sec. 3, and examples in Sec. 4.
This implies that for a vector 
x=(alpha_1,alpha_2,...,alpha_m,...,alpha_N)^T with geometrical parameters 
alpha_m, and unitary vectors e_m, m=1,2,N, I need to perform the 
assembly for f(x+eps*e_m), f(x).
I hope this clarifies my aim with this email.

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.

On Sunday, 19 January 2020 23:50:37 UTC+1, Jean-Paul Pelteret wrote:
>
> 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
>
> On 16 Jan 2020, at 13:19, Juan Carlos Araujo Cabarcas  > 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 

Re: [deal.II] Different shape representations with manifolds on the same triangulation

[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

> On 16 Jan 2020, at 13:19, Juan Carlos Araujo Cabarcas  
> 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 
>  
> 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 

Re: [deal.II] Different shape representations with manifolds on the same triangulation

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


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 
>  
> 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
>
> On 15 Jan 2020, at 14:46, Juan Carlos Araujo Cabarcas  > 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 _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 polar;
> TransfiniteInterpolationManifold inner_manifold;
>
> Triangulation triangulation;
>
> hp::DoFHandlerdof_handler;
> hp::FECollectionfe_collection;
> hp::MappingCollection mapping_collection;
>
> PETScWrappers::SparseMatrixsystem_matrix;
> PETScWrappers::MPI::Vector solution;
>
> // Constructor
> manifold ( alpha), dof_handler (triangulation),...
>
> // setup_system
> typename hp::DoFHandler::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 ( 

[deal.II] 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 _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 polar;
TransfiniteInterpolationManifold inner_manifold;

Triangulation triangulation;

hp::DoFHandlerdof_handler;
hp::FECollectionfe_collection;
hp::MappingCollection mapping_collection;

PETScWrappers::SparseMatrixsystem_matrix;
PETScWrappers::MPI::Vector solution;

// Constructor
manifold ( alpha), dof_handler (triangulation),...

// setup_system
typename hp::DoFHandler::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 hp_fe_v ( mapping_collection,

fe_collection,
quadrature_collection,
update_values|  update_gradients |
update_quadrature_points  | 
 update_JxW_values);

FullMatrix   cell_system;
Vectorcell_rhs;

std::vector local_dof_indices;

typename hp::DoFHandler::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 _v = hp_fe_v.get_present_fe_values ();

for (unsigned int q_point=0; q_point x = fe_v.quadrature_point (q_point);
  
  for (unsigned int i=0; ihttp://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/92fc053d-ca4e-4830-bca0-2f211a9ba384%40googlegroups.com.