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

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

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

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

[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 &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 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 &fe_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.

Re: [deal.II] Re: optimization, SolverBFGS

[Juan Carlos Araujo Cabarcas]

Thanks for the help, now it is running nicely!

On Monday, 28 October 2019 15:05:19 UTC+1, Bruno Turcksin wrote:
>
> There is a small bug in the lambda: 
>
> Le lun. 28 oct. 2019 à 09:36, Juan Carlos Araujo Cabarcas 
> > a écrit : 
> >   auto my_lambda = [&polar](Vector a, Vectorb) { 
> >   return polar.objective_fun(a,b); 
> >   };  //polar.Callback (a,b);}; 
>
> You want the vector to be passed by reference not by value, otherwise 
> the result cannot be updated. So it should be: 
> auto my_lambda = [&polar](const Vector &a, Vector &b) { 
>return polar.objective_fun(a,b); 
> }; 
>
> Bruno 
>

-- 
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/110bfff8-6341-42df-a010-1c1a86e50435%40googlegroups.com.

Re: [deal.II] Re: optimization, SolverBFGS

[Juan Carlos Araujo Cabarcas]

Hi, thanks! that worked by using polar.objective_fun(a,b), instead of 
Callback.
Now, it seems that SolverBFGS is not starting, and terminates in the first 
function call! 
However, the optimization should not be done since |grad x| is not small. 

This is what I do:

 Vector X (N), Dx (N);
  
  SolverControlresidual_control (100, 1e-10, true, 
true);
  residual_control.enable_history_data();
  
  SolverBFGS > opt (residual_control, 
  
SolverBFGS >::AdditionalData(50, true) );
  
  Shape_compute  polar(par);
  
  auto my_lambda = [&polar](Vector a, Vectorb) {
  return polar.objective_fun(a,b);
  };  //polar.Callback (a,b);};
opt.solve ( my_lambda, X );
  
  std::cout << "  converged in " << residual_control.last_step()
  << " iterations" << std::endl;
  
  polar.objective_fun(X,Dx);  
  printf("\nThe norm of Dx is %.5f", Dx.l1_norm() );

and I get 
  converged in 0 iterations
with The norm of Dx is 0.85597


On Monday, October 28, 2019 at 1:24:22 PM UTC+1, Bruno Turcksin wrote:
>
> Here is a simple example: 
> https://wandbox.org/permlink/fbcldo9PIiHl6tfI  In your case, try 
> something like 
>
> auto my_lambda = [&polar](Vector a, Vectorb) {return 
> polar.Callback(a,b);}; 
> opt.solve ( my_lambda, X ); 
>
> Bruno 
>
> Le lun. 28 oct. 2019 à 06:30, Juan Carlos Araujo Cabarcas 
> > a écrit : 
> > 
> > Thanks Bruno for your prompt reply, 
> > 
> > I got the important part of the code compiling, however since I am new 
> to this lambda (functional) functions, I do not know how to pass the new 
> definitions to SolverBFGS. 
> > 
> > Looking at your suggestion I do the following: 
> > 
> > class Shape_compute { 
> > public: 
> > Parameters par; 
> > 
> > Shape_compute (Parameters epar):par(epar) { 
> > Register([=]( Vector ex, Vector dx ){ return 
> objective_fun (ex,dx); }); 
> > } 
> > 
> > void Register(std::function, Vector)> 
> Callback) {} //fun = objective_fun; 
> > 
> > double objective_fun ( Vector ex, Vector dx) { 
> > 
> > vector x (par.design_N), grad (par.design_N); 
> > for (unsigned int i=0; i < x.size (); i++) { 
> > x [i] = ex(i); 
> > } 
> > Adaptive::DtN_Helmholtz<2> fem (par, x, true); 
> > fem.run (); 
> > double objective_function = fem.objective_function; 
> > 
> > for (unsigned int i=0; i < par.design_N; i++) { 
> > if (i < par.design_N) { 
> > vector aux = x; aux [i] += h; 
> > Adaptive::DtN_Helmholtz<2> fem (par, aux, false); 
> > fem.run (); 
> > grad [i] = (fem.objective_function - 
> objective_function )/h; 
> > } 
> > dx(i) = grad [i]; 
> > } 
> > iteration++; 
> > return objective_function; 
> > } 
> > }; 
> > 
> > 
> > Then, when calling the optimization routine: 
> > 
> >   Vector X (N); 
> > 
> >   SolverControlresidual_control (N, 1e-7); 
> >   //SolverBFGS::AdditionalData   data (10, true); 
> > 
> >   SolverBFGS > opt (residual_control);//, data); 
> > 
> >   Shape_compute  polar(par); 
> >   opt.solve ( polar.Callback, X ); // If this line is commented, 
> it compiles fine! 
> > 
> > 
> > I would appreciate guidance on how to achieve this, 
> > 
> > Thanks in advance, 
> >   Juan Carlos Araújo 
>

-- 
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/349c3304-a316-4364-a230-a5b847f7d903%40googlegroups.com.

[deal.II] Re: optimization, SolverBFGS

[Juan Carlos Araujo Cabarcas]

Thanks Bruno for your prompt reply,

I got the important part of the code compiling, however since I am new to 
this lambda (functional) functions, I do not know how to pass the new 
definitions to SolverBFGS.

Looking at your suggestion I do the following:

class Shape_compute {
public:
Parameters par;

Shape_compute (Parameters epar):par(epar) {
Register([=]( Vector ex, Vector dx ){ return 
objective_fun (ex,dx); });
}

void Register(std::function, Vector)> 
Callback) {} //fun = objective_fun;

double objective_fun ( Vector ex, Vector dx) {

vector x (par.design_N), grad (par.design_N);
for (unsigned int i=0; i < x.size (); i++) {
x [i] = ex(i);
}
Adaptive::DtN_Helmholtz<2> fem (par, x, true);
fem.run ();
double objective_function = fem.objective_function;

for (unsigned int i=0; i < par.design_N; i++) {
if (i < par.design_N) {
vector aux = x; aux [i] += h;
Adaptive::DtN_Helmholtz<2> fem (par, aux, false);
fem.run (); 
grad [i] = (fem.objective_function - objective_function 
)/h;
}
dx(i) = grad [i];
}
iteration++;
return objective_function;
}
}; 


Then, when calling the optimization routine:

  Vector X (N);
  
  SolverControlresidual_control (N, 1e-7);
  //SolverBFGS::AdditionalData   data (10, true);
  
  SolverBFGS > opt (residual_control);//, data);
  
  Shape_compute  polar(par);
  opt.solve ( polar.Callback, X ); // If this line is commented, it 
compiles fine!


I would appreciate guidance on how to achieve this,

Thanks in advance,
  Juan Carlos Araújo




On Friday, 25 October 2019 17:55:03 UTC+2, Bruno Turcksin wrote:
>
> Juan,
>
> Take a look at 
> https://stackoverflow.com/questions/17131768/how-to-directly-bind-a-member-function-to-an-stdfunction-in-visual-studio-11
>  
> I advise using a lambda like the second reply suggests
>
> Best,
>
> Bruno
>
> On Friday, October 25, 2019 at 10:53:08 AM UTC-4, Juan Carlos Araujo 
> Cabarcas wrote:
>>
>> Dear all,
>>
>> I would like to use SolverBFGS from deal.II/optimization/ in my project, 
>> and I try to follow the documentation in
>> https://www.dealii.org/current/doxygen/deal.II/classSolverBFGS.html
>>
>> First, I would like to mention that the documentation would really 
>> improve by adding a minimal example on how to use SolverBFGS.
>>
>> In my project, I cannot just define an independent function to be pass to 
>> SolverBFGS.
>> Instead, I define a class ShapePolar, which needs to be initialized with 
>> some mesh parameters par, as in
>>
>> ShapePolar polar (par);
>> 
>> then I have implemented the function 
>>
>> double ShapePolar::fun (const Vector, Vector)
>>
>> that computes a scalar value of an objective function and its gradient 
>> Vector, when we pass X.
>>
>>
>> Then, my attempt, after looking at the documentation, is to do something 
>> like this:
>>
>>   Vector X (N);
>>   
>>   SolverControlresidual_control (N, 1e-7);
>>   //SolverBFGS::AdditionalData   data (10, true);
>>
>>   SolverBFGS > opt (residual_control);   //, data);
>>
>>   ShapePolar  polar (par);
>>   std::function, Vector)>  fun = 
>> polar.fun;
>>
>>   opt.solve ( fun, X );
>>
>> After several failed attempts (syntax errors), I realize that I am 
>> missing something fundamental here and could use some guidance.
>> I would appreciate any hint on how to achieve using SolverBFGS for the 
>> setting described above.
>>
>> Kind regards, 
>> Juan C. Araújo Cabarcas
>>
>>

-- 
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/599f60b4-a71e-4d6f-9b4c-f89c48078dc4%40googlegroups.com.

[deal.II] optimization, SolverBFGS

[Juan Carlos Araujo Cabarcas]

Dear all,

I would like to use SolverBFGS from deal.II/optimization/ in my project, 
and I try to follow the documentation in
https://www.dealii.org/current/doxygen/deal.II/classSolverBFGS.html

First, I would like to mention that the documentation would really improve 
by adding a minimal example on how to use SolverBFGS.

In my project, I cannot just define an independent function to be pass to 
SolverBFGS.
Instead, I define a class ShapePolar, which needs to be initialized with 
some mesh parameters par, as in

ShapePolar polar (par);

then I have implemented the function 

double ShapePolar::fun (const Vector, Vector)

that computes a scalar value of an objective function and its gradient 
Vector, when we pass X.


Then, my attempt, after looking at the documentation, is to do something 
like this:

  Vector X (N);
  
  SolverControlresidual_control (N, 1e-7);
  //SolverBFGS::AdditionalData   data (10, true);

  SolverBFGS > opt (residual_control);   //, data);

  ShapePolar  polar (par);
  std::function, Vector)>  fun = 
polar.fun;

  opt.solve ( fun, X );

After several failed attempts (syntax errors), I realize that I am missing 
something fundamental here and could use some guidance.
I would appreciate any hint on how to achieve using SolverBFGS for the 
setting described above.

Kind regards, 
Juan C. Araújo Cabarcas

-- 
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/6b633322-dd52-4418-933d-86a798265936%40googlegroups.com.

Re: [deal.II] How to manually create sparsity pattern for PETSc sparsity matrix in parallel

[Juan Carlos Araujo Cabarcas]

Hey Pai Liu,

At some point in my research I had to do something similar to what you want 
to do. In the application of the so-called Dirichlet-to-Neumann maps for 
Helmholtz problems, one has to assemble a dense block corresponding to the 
boundary degrees of freedom, so I had to include extra entries in my FE 
matrices. This is what I did (omitting some non-related details):

PETScWrappers::SparseMatrix   Q;

std::vector map_dofs_boundary;
DoFTools::map_dof_to_boundary_indices (dof_handler,

map_dofs_boundary);

std::vector dofs_boundary_lenghts ( map_dofs_boundary.size() 
);

Q.reinit (par.dofs, par.dofs, dofs_boundary_lenghts);
MatSetOption(Q, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

FE assembly routine ...

in a loop i, j ...
Q.set ( boundary_idx[i], boundary_idx[j], my_value (i,j) );

Q.compress (VectorOperation::add);

and that works for me.

I hope this may be of help to you.

Juan Carlos Araújo Cabarcas
Doctor candidate in Computational Science with emphasis in mathematics,
Umeå Universitet


El miércoles, 3 de abril de 2019, 0:47:09 (UTC+2), Pai Liu escribió:
>
> Hi Wolfgang,
>  
>
>> I would start by asking where the nonzero entries are going to lie. What 
>> kind of operator do you plan on implementing for which you don't know 
>> where the nonzeros are? 
>>
>
> I want to multiply this matrix with a PETSc vector (which is also stored 
> in parallel).
> When I initialize the PETSc sparsity matrix, I don't know where the 
> nonzeros are.
> Then I loop over the locally owned cells on each processor, and guadually 
> know for each row of the matrix the locations and values of the nonzeros. 
> Do you mean that I can't directly use the PETSc_sparsity_matrix.add() 
> function to insert nonzero values into the matrix without telling the 
> sparsity pattern the nonzeros locations?
>
> *If I have to tell the sparsity pattern where the nonzeros are manually, 
> what are the necessary procedures and functions I should follow and use 
> (since the matrix is not a system matrix, I can't use the 
> DoFTools::make_sparsity_pattern function with an dofhandler as input)?*
>
> *For example, for row 0, how can I tell the sparsity pattern the elements 
> at (0,3) and (0,5) arenonzeros? I search the documentation, and find maybe 
> I should use dynamic_spasity_pattern? I am so confused.*
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: shape_value, shape_grad in arbitrary points within quad that are not quadrature points

[Juan Carlos Araujo Cabarcas]

Thanks a lot for the pointers! 

El viernes, 17 de agosto de 2018, 22:57:24 (UTC+2), Wolfgang Bangerth 
escribió:
>
> On 08/17/2018 02:56 PM, Juan Carlos Araujo Cabarcas wrote: 
> > Hej, thanks for your prompt answer! 
> > Yes, your interpretation of my poor description is correct. So then I 
> just 
> > need to create a new FEValues within the cell loop containing the points 
> where 
> > I want to evaluate the shape functions, right? How do I add this points? 
> > 
> > Also if the points are in the physical quad, do I need to map them to 
> the 
> > reference quad [-1,1]x[-1,1] before adding them, or I just add the 
> points 
> > directly to the FEValues? 
>
> Yes. 
>
> I think you probably want to read through the implementation of the 
> VectorTools::point_value() function, which does exactly what you want. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: shape_value, shape_grad in arbitrary points within quad that are not quadrature points

[Juan Carlos Araujo Cabarcas]

Hej, thanks for your prompt answer! 
Yes, your interpretation of my poor description is correct. So then I just 
need to create a new FEValues within the cell loop containing the points 
where I want to evaluate the shape functions, right? How do I add this 
points? 

Also if the points are in the physical quad, do I need to map them to the 
reference quad [-1,1]x[-1,1] before adding them, or I just add the points 
directly to the FEValues?

El viernes, 17 de agosto de 2018, 18:18:41 (UTC+2), Daniel Arndt escribió:
>
> Juan,
>
> I would like to ask what would be the best/fastest way to evaluate 
> *fe_v.shape_value 
>> (i,x), **fe_v.shape_grad* *(i,x)*, in points that are not quadrature 
>> points but that are guaranteed to lie within the curved cell/quad.
>>
> If I understand the situation correctly, you have a FEValues object that 
> is initialized on a cell and now you want to evaluate it at an additional 
> point that is not included in the mapped quadrature points.
> Is that the case? If so, I would either add that point to the existing 
> FEValues object or build a separate one that includes the quadrature point 
> (with respect to the unit cell) you want to evaluate the shape functions at.
> It is sufficient to initialize this separate object for cells for which 
> you actually need it.
> Maybe, I misunderstand what you want. Can you specify the issue you are 
> trying to solve?
>
> Best,
> Daniel
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: installation error

[Juan Carlos Araujo Cabarcas]

Hi, when I install PETSc I allow the configuration process to download the 
necessary libraries it needs: 

./config/configure.py --with-shared=1 --with-x=0 --with-mpi=1 
--with-scalar-type=complex --download-mumps --download-metis 
--download-parmetis --download-superlu_dist --download-blacs 
--download-scalapack --download-superlu_dist

>From the dealii detailed.log I see that:

MPI_VERSION = 2.1
OMPI_VERSION = 1.6.5

>From the terminal:

$ mpicc --version
gcc (Ubuntu 4.8.4-2ubuntu1~14.04.3) 4.8.4

So it seems I need to upgrade to MPI 3.0.
Thanks, for having a look at this.

El jueves, 25 de enero de 2018, 12:29:53 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 01/24/2018 07:34 AM, Juan Carlos Araujo Cabarcas wrote: 
> > Please find the file: detailed.log attached. 
> > 
> > El martes, 23 de enero de 2018, 17:02:14 (UTC-5), Wolfgang Bangerth 
> > escribió: 
> > 
> > On 01/23/2018 02:13 PM, Bruno Turcksin wrote: 
> >  > 
> >  > mypath/dealii/source/lac/scalapack.cc:243:91: error: there 
> > are no 
> >  > arguments to ‘MPI_Comm_create_group’ that depend on a 
> > template parameter, 
> >  > so a declaration of ‘MPI_Comm_create_group’ must be available 
> > [-fpermissive] 
> >  > ierr = MPI_Comm_create_group(MPI_COMM_WORLD, 
> > group_union, 5, 
> >  > &mpi_communicator_union); 
>
> It confused me that you are compiling the scalapack file but get an 
> error message that a particular MPI function was not found. This would 
> ordinarily suggest that either scalapack.cc is missing an #include of 
>  (which is not the case) or that your installation is not 
> configured with MPI (which is the case for you). So the error did not 
> make sense to me at first. 
>
> But it turns out that the call in question, MPI_Comm_create_group is a 
> function that was only added to MPI in version 3.0. Apparently all of us 
> use an MPI installation that is sufficiently up to date. What is the 
> version you use? 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: installation error

[Juan Carlos Araujo Cabarcas]

Hi, when I install PETSc I allow the configuration process to download the 
necessary libraries it needs: 

./config/configure.py --with-shared=1 --with-x=0 --with-mpi=1 
--with-scalar-type=complex --download-mumps --download-metis 
--download-parmetis --download-superlu_dist --download-blacs 
--download-scalapack --download-superlu_dist

>From the dealii detailed.log I see that:

MPI_VERSION = 2.1
OMPI_VERSION = 1.6.5

>From the terminal:

$ mpicc --version
gcc (Ubuntu 4.8.4-2ubuntu1~14.04.3) 4.8.4

Shall I update my system libraries?

El jueves, 25 de enero de 2018, 12:29:53 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 01/24/2018 07:34 AM, Juan Carlos Araujo Cabarcas wrote: 
> > Please find the file: detailed.log attached. 
> > 
> > El martes, 23 de enero de 2018, 17:02:14 (UTC-5), Wolfgang Bangerth 
> > escribió: 
> > 
> > On 01/23/2018 02:13 PM, Bruno Turcksin wrote: 
> >  > 
> >  > mypath/dealii/source/lac/scalapack.cc:243:91: error: there 
> > are no 
> >  > arguments to ‘MPI_Comm_create_group’ that depend on a 
> > template parameter, 
> >  > so a declaration of ‘MPI_Comm_create_group’ must be available 
> > [-fpermissive] 
> >  > ierr = MPI_Comm_create_group(MPI_COMM_WORLD, 
> > group_union, 5, 
> >  > &mpi_communicator_union); 
>
> It confused me that you are compiling the scalapack file but get an 
> error message that a particular MPI function was not found. This would 
> ordinarily suggest that either scalapack.cc is missing an #include of 
>  (which is not the case) or that your installation is not 
> configured with MPI (which is the case for you). So the error did not 
> make sense to me at first. 
>
> But it turns out that the call in question, MPI_Comm_create_group is a 
> function that was only added to MPI in version 3.0. Apparently all of us 
> use an MPI installation that is sufficiently up to date. What is the 
> version you use? 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: installation error

[Juan Carlos Araujo Cabarcas]

Please find the file: detailed.log attached.

El martes, 23 de enero de 2018, 17:02:14 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 01/23/2018 02:13 PM, Bruno Turcksin wrote: 
> > 
> > mypath/dealii/source/lac/scalapack.cc:243:91: error: there are no 
> > arguments to ‘MPI_Comm_create_group’ that depend on a template 
> parameter, 
> > so a declaration of ‘MPI_Comm_create_group’ must be available 
> [-fpermissive] 
> > ierr = MPI_Comm_create_group(MPI_COMM_WORLD, group_union, 5, 
> > &mpi_communicator_union); 
> > 
> > Do you need scalapack? If you don't, just add DEAL_II_WITH_SCALAPACK=OFF 
>
> Out of curiosity, can you also send the output of the detailed.log file in 
> your build directory? 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
###
#
#  deal.II configuration:
#CMAKE_BUILD_TYPE:   DebugRelease
#BUILD_SHARED_LIBS:  ON
#CMAKE_INSTALL_PREFIX:   /home/ju4nk4/Soft/dealii
#CMAKE_SOURCE_DIR:   /home/ju4nk4/Soft/dealii
#(version 9.0.0-pre, shortrev 55c9e45)
#CMAKE_BINARY_DIR:   /home/ju4nk4/Soft/dealii/build
#CMAKE_CXX_COMPILER: GNU 4.8.4 on platform Linux x86_64
#/usr/bin/c++
#CMAKE_C_COMPILER:   /usr/bin/cc
#CMAKE_Fortran_COMPILER: /usr/bin/gfortran
#CMAKE_GENERATOR:Unix Makefiles
#
#  Base configuration (prior to feature configuration):
#DEAL_II_CXX_FLAGS:-pedantic -fPIC -Wall -Wextra 
-Wpointer-arith -Wwrite-strings -Wsynth -Wsign-compare -Wswitch 
-Woverloaded-virtual -Wno-deprecated-declarations -Wno-literal-suffix -std=c++11
#DEAL_II_CXX_FLAGS_RELEASE:-O2 -funroll-loops -funroll-all-loops 
-fstrict-aliasing -Wno-unused-local-typedefs
#DEAL_II_CXX_FLAGS_DEBUG:  -Og -ggdb -Wa,--compress-debug-sections
#DEAL_II_LINKER_FLAGS: -Wl,--as-needed -rdynamic -fuse-ld=gold
#DEAL_II_LINKER_FLAGS_RELEASE: 
#DEAL_II_LINKER_FLAGS_DEBUG:   -ggdb
#DEAL_II_DEFINITIONS:  
#DEAL_II_DEFINITIONS_RELEASE:  
#DEAL_II_DEFINITIONS_DEBUG:DEBUG
#DEAL_II_USER_DEFINITIONS: 
#DEAL_II_USER_DEFINITIONS_REL: 
#DEAL_II_USER_DEFINITIONS_DEB: DEBUG
#DEAL_II_INCLUDE_DIRS  
#DEAL_II_USER_INCLUDE_DIRS:
#DEAL_II_BUNDLED_INCLUDE_DIRS: 
#DEAL_II_LIBRARIES:m
#DEAL_II_LIBRARIES_RELEASE:
#DEAL_II_LIBRARIES_DEBUG:  
#
#  Configured Features (DEAL_II_ALLOW_BUNDLED = ON, DEAL_II_ALLOW_AUTODETECTION 
= ON):
#  ( DEAL_II_WITH_64BIT_INDICES = OFF )
#  ( DEAL_II_WITH_ADOLC = OFF )
#DEAL_II_WITH_ARPACK set up with external dependencies
#ARPACK_LINKER_FLAGS = 
#ARPACK_LIBRARIES = 
/usr/lib/libarpack.so;/usr/lib/liblapack.so;/usr/lib/libblas.so;gfortran;quadmath;m;c;/usr/lib/openmpi/lib/libmpi.so;dl;/usr/lib/x86_64-linux-gnu/libhwloc.so
#  ( DEAL_II_WITH_ASSIMP = OFF )
#DEAL_II_WITH_BOOST set up with bundled packages
#BOOST_CXX_FLAGS = -Wno-unused-local-typedefs
#BOOST_BUNDLED_INCLUDE_DIRS = 
/home/ju4nk4/Soft/dealii/bundled/boost-1.62.0/include
#BOOST_LIBRARIES = rt
#  ( DEAL_II_WITH_CUDA = OFF )
#  ( DEAL_II_WITH_CXX14 = OFF )
#  ( DEAL_II_WITH_CXX17 = OFF )
#DEAL_II_WITH_GMSH set up with external dependencies
#GMSH_EXE = /usr/bin/gmsh
#  ( DEAL_II_WITH_GSL = OFF )
#  ( DEAL_II_WITH_HDF5 = OFF )
#DEAL_II_WITH_LAPACK set up with external dependencies
#LAPACK_WITH_64BIT_BLAS_INDICES = OFF
#LAPACK_LINKER_FLAGS = 
#LAPACK_LIBRARIES = 
/usr/lib/liblapack.so;/usr/lib/libblas.so;gfortran;quadmath;m;c
#DEAL_II_WITH_METIS set up with external dependencies
#METIS_VERSION = 5.1.0
#METIS_DIR = /home/ju4nk4/Soft/metis-5.1.0
#METIS_INCLUDE_DIRS = /home/ju4nk4/Soft/metis-5.1.0/include
#METIS_USER_INCLUDE_DIRS = /home/ju4nk4/Soft/metis-5.1.0/include
#METIS_LIBRARIES = 
/home/ju4nk4/Soft/metis-5.1.0/build/Linux-x86_64/libmetis/libmetis.a;/usr/lib/openmpi/lib/libmpi.so;dl;/usr/lib/x86_64-linux-gnu/libhwloc.so
#DEAL_II_WITH_MPI set up with external dependencies
#MPI_V

[deal.II] installation error

[Juan Carlos Araujo Cabarcas]

Dear all,

I am trying to install deal.II from the GIT repository with the following 
features:

petsc_ver='3.6.0';
trilinos_ver='12.4.2';

git clone https://github.com/dealii/dealii.git dealii

cmake \
-DTRILINOS_DIR=${install_dir}/trilinos-${trilinos_ver} \
-DP4EST_DIR=${install_dir}/FAST \
-DDEAL_II_WITH_METIS=ON \
-DMETIS_DIR=$METIS_DIR \
-DDEAL_II_WITH_MPI=ON \
-DDEAL_II_WITH_THREADS=OFF \
-DDEAL_II_WITH_UMFPACK=ON \
-DDEAL_II_WITH_LAPACK=ON \
-DDEAL_II_WITH_PETSC=ON \
-DPETSC_ARCH=$PETSC_ARCH \
-DPETSC_DIR=$PETSC_DIR \
-DDEAL_II_WITH_SLEPC=ON \
-DSLEPC_DIR=$SLEPC_DIR \
-DDEAL_II_WITH_P4EST=ON \
-DDEAL_II_WITH_ARPACK=ON \
-DDEAL_II_WITH_TRILINOS=ON \
-DCMAKE_INSTALL_PREFIX=${install_dir}/dealii ${install_dir}/dealii;

make -j${virtual_proc} install;

I get the following installation error:

[ 50%] Building CXX object 
source/fe/CMakeFiles/obj_fe_debug.dir/fe_poly.cc.o
mypath/dealii/source/lac/scalapack.cc: In member function ‘void 
dealii::ScaLAPACKMatrix::copy_to(dealii::ScaLAPACKMatrix&)
 
const’:
mypath/dealii/source/lac/scalapack.cc:243:91: error: there are no arguments 
to ‘MPI_Comm_create_group’ that depend on a template parameter, so a 
declaration of ‘MPI_Comm_create_group’ must be available [-fpermissive]
   ierr = MPI_Comm_create_group(MPI_COMM_WORLD, group_union, 5, 
&mpi_communicator_union);

   
^
mypath/dealii/source/lac/scalapack.cc:243:91: note: (if you use 
‘-fpermissive’, G++ will accept your code, but allowing the use of an 
undeclared name is deprecated)
mypath/dealii/source/lac/scalapack.cc: In instantiation of ‘void 
dealii::ScaLAPACKMatrix::copy_to(dealii::ScaLAPACKMatrix&)
 
const [with NumberType = double]’:
mypath/dealii/source/lac/scalapack.cc:760:16:   required from here
mypath/dealii/source/lac/scalapack.cc:243:91: error: 
‘MPI_Comm_create_group’ was not declared in this scope
mypath/dealii/source/lac/scalapack.cc: In instantiation of ‘void 
dealii::ScaLAPACKMatrix::copy_to(dealii::ScaLAPACKMatrix&)
 
const [with NumberType = float]’:
mypath/dealii/source/lac/scalapack.cc:761:16:   required from here
mypath/dealii/source/lac/scalapack.cc:243:91: error: 
‘MPI_Comm_create_group’ was not declared in this scope
[ 50%] Building CXX object 
source/grid/CMakeFiles/obj_grid_debug.dir/grid_reordering.cc.o
make[2]: *** [source/lac/CMakeFiles/obj_lac_debug.dir/scalapack.cc.o] Error 
1
make[1]: *** [source/lac/CMakeFiles/obj_lac_debug.dir/all] Error 2
make[1]: *** Waiting for unfinished jobs

Last time I had errors was due to Trilinos version, but I updated that one 
to the current in use.
Any ideas for resolving the error?

Aditionally, how can I set from GIT the last stable (tested) version?

Thanks in advance,
Juan Carlos Araújo, 
Umeå Universitet

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Error with Transfinite Interpolation

[Juan Carlos Araujo Cabarcas]

Thanks a lot for looking at this issue so quickly!

No problem, you can use it. That was the idea of sharing the minimal 
example! I am always happy to contribute with the library.

So, does this mean that the way I colorize is correct? if there is a 
better/simpler way plz let me know.
Is the solution merged with the latest master branch? I would like to have 
a working library for this particular example.

El jueves, 18 de enero de 2018, 18:39:51 (UTC-5), David Wells escribió:
>
> I figured it out: we have an optimization that assumes in one particular 
> place that, if we have eight points, they must be the vertices of a cube. 
> However, in 2D, we use eight points to place new points on quadrilaterals: 
> assuming these are the vertices of a cube is wrong. The original assumption 
> created a bad initial value for a quasi-newton solver, which then failed to 
> converge.
>
> I got rid of the cube assumption. see 
> https://github.com/dealii/dealii/pull/5761 for more information.
>
> Do you mind if I use a modified version of your code for a new test in 
> deal.II?
>
> Thanks,
> David Wells
>
> On Wednesday, January 17, 2018 at 3:05:48 PM UTC-5, Juan Carlos Araujo 
> Cabarcas wrote:
>>
>> Dear all,
>>
>> Since the inclusion of Transfinite interpolation, I have been successful 
>> on working with this powerful technique in my research. I had coded a mesh 
>> implementing concentric circles, where the inner most is shifted a small 
>> distance s. All concentric circles are labeled 100+i, where i is in a loop 
>> on all circles.
>>
>> The mesh was working after 4/Apr/17 when I added the following entry in 
>> the forum:
>> https://groups.google.com/forum/#!topic/dealii/hCZqv9g6mKk
>>
>> I installed the development version of dealii on the 17/nov/17. The 
>> details of the installation are also in the forum:
>> https://groups.google.com/forum/#!topic/dealii/ee2w2X987ig
>>
>> and recently I went back and tried to run thhis code, obtaining the error 
>> at the end of the email.
>>
>> I prepared a minimal example that includes the definition of my grid. It 
>> includes how I colorize each face as explained before, and the rest are 
>> colorized "1" as the documentation suggests.
>>
>> The symptoms are the following:
>> - If I use refinement=0, everything works and the mesh looks quite nice 
>> in zeros.vtk and surface with lines mode in Paraview.
>> - For any higher refinement, it gives the error below.
>>
>> Maybe there has been a major change in TransfiniteInterpolation since 
>> Apr/2017, or maybe the way I colorize is not good anymore for some weird 
>> reason. Any ideas or comments would be greatly appreciated.
>>
>> Juan Carlos Araújo,
>> Umeå Universitet.
>>
>>
>>
>> terminate called after throwing an instance of 'dealii::Mapping<2, 
>> 2>::ExcTransformationFailed'
>>   what():  
>> 
>> An error occurred in line <1554> of file 
>>  in function
>> typename dealii::Triangulation::cell_iterator 
>> dealii::TransfiniteInterpolationManifold> spacedim>::compute_chart_points(const dealii::ArrayView> dealii::Point >&, dealii::ArrayView >) const 
>> [with int dim = 2; int spacedim = 2; typename dealii::Triangulation> spacedim>::cell_iterator = dealii::TriaIterator 
>> >]
>> The violated condition was: 
>> false
>> Additional information: 
>> 
>> 
>>
>> Aborted (core dumped)
>>
>>
>>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Error with Transfinite Interpolation

[Juan Carlos Araujo Cabarcas]

Dear all,

Since the inclusion of Transfinite interpolation, I have been successful on 
working with this powerful technique in my research. I had coded a mesh 
implementing concentric circles, where the inner most is shifted a small 
distance s. All concentric circles are labeled 100+i, where i is in a loop 
on all circles.

The mesh was working after 4/Apr/17 when I added the following entry in the 
forum:
https://groups.google.com/forum/#!topic/dealii/hCZqv9g6mKk

I installed the development version of dealii on the 17/nov/17. The details 
of the installation are also in the forum:
https://groups.google.com/forum/#!topic/dealii/ee2w2X987ig

and recently I went back and tried to run thhis code, obtaining the error 
at the end of the email.

I prepared a minimal example that includes the definition of my grid. It 
includes how I colorize each face as explained before, and the rest are 
colorized "1" as the documentation suggests.

The symptoms are the following:
- If I use refinement=0, everything works and the mesh looks quite nice in 
zeros.vtk and surface with lines mode in Paraview.
- For any higher refinement, it gives the error below.

Maybe there has been a major change in TransfiniteInterpolation since 
Apr/2017, or maybe the way I colorize is not good anymore for some weird 
reason. Any ideas or comments would be greatly appreciated.

Juan Carlos Araújo,
Umeå Universitet.



terminate called after throwing an instance of 'dealii::Mapping<2, 
2>::ExcTransformationFailed'
  what():  

An error occurred in line <1554> of file 
 in function
typename dealii::Triangulation::cell_iterator 
dealii::TransfiniteInterpolationManifold::compute_chart_points(const dealii::ArrayView >&, dealii::ArrayView >) const 
[with int dim = 2; int spacedim = 2; typename dealii::Triangulation::cell_iterator = dealii::TriaIterator 
>]
The violated condition was: 
false
Additional information: 



Aborted (core dumped)


-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 1999 - 2013 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -
 *
*/

#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 

// #include "../../grid/arbitrary_resonators.h"

#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 


using namespace dealii;
using namespace std;

// GRID DEFINITION -

struct Geom_parameters {
  
  unsigned int n_scatterers, scatterers_coarse_cells, 
  	n_balls, coated_coarse_cells, defective_index;
  	
  unsigned int manifold_label,layers_label,bc_label;
  
  double coat_radius;
  bool coated=false;
  bool defect=false;
  bool defect_n2=false;
  vector< Point<2> > ball_centers;
  vector radius;
  vector radius_PML;
  vector ball_labels;
  
  string name;
};

void concentric_disks (	Triangulation<2> &tria, double s, vector x,
Geom_parameters &gp, bool mesh_out)
{
	double r = x[0], d = 0.5*x[0], q=1.0/sqrt(2.0); //l = x[1], q: corner points factor
	
	const unsigned int nlay = x.size()-1, vlayer=8, lcells=8, 
			n_vert = (2+nlay)*vlayer+1, n_cell = 3+(1+nlay)*lcells+1; // 19
  
  
  bool info = false;	//, plot_eps = true, for printing cells, faces and vertex values
	
	int fc = 4; //fv = 1,  f: fixed, displacement of index on the vertexes and cells ... inner objects
  
	// geometric information---
			gp.n_scatterers = 1;
			gp.scatterers_coarse_cells = 12;
			gp.n_balls = x.size(); // the resonator and all other centered at the origin
			gp.ball_centers.resize( gp.n_balls );
			
			gp.layers_label=1; 
			gp.manifold_label=10;
			gp.bc_label=0;
			
			gp.radius.resize(gp.n_balls);
			
			gp.radius[0]=x[0];
			gp.ball_centers[0] = Point<2>( s, 0.0 );
			
			for (unsigned int k = 1; k < x.size(); k++) {
gp.radius[k]=x[k];
gp.ball_centers[k] =

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Ahh I see. Thanks to all of you for all the help and support on this thread!

El sábado, 9 de diciembre de 2017, 8:42:38 (UTC-5), Daniel Arndt escribió:
>
> Juan,
>
> take a look at the documentation 
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceMatrixTools.html#a6d63cb21abc6c954557c603478ba5f7d>
>  of 
> that function. It tells you that the last parameter needs to be "false", 
> i.e.
> this function is only implemented for the case in which you don't 
> eliminate the columns corresponding to boundary DoFs.
>
> Best,
> Daniel
>
> Am Samstag, 9. Dezember 2017 00:48:04 UTC+1 schrieb Juan Carlos Araujo 
> Cabarcas:
>>
>> Thanks again, 
>>
>> Now it passes the *interpolate_boundary_values *call. However, it is 
>> conflicting now with:
>>
>>   MatrixTools::apply_boundary_values (bval,
>>   
>>   system_matrix,
>>   
>>   solution,
>>   
>>   system_rhs);
>>
>> To be precise I write the output below. 
>>
>> Thanks in advance,
>> Juan Carlos Araújo, Umeå University
>>
>>
>> 
>> An error occurred in line <70> of file 
>>  in function
>> void dealii::MatrixTools::apply_boundary_values(const 
>> std::map >&, 
>> dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::VectorBase&, 
>> dealii::PETScWrappers::VectorBase&, bool)
>> The violated condition was: 
>> eliminate_columns == false
>> Additional information: 
>> You are trying to use functionality in deal.II that is currently not 
>> implemented. In many cases, this indicates that there simply didn't appear 
>> much of a need for it, or that the author of the original code did not have 
>> the time to implement a particular case. If you hit this exception, it is 
>> therefore worth the time to look into the code to find out whether you may 
>> be able to implement the missing functionality. If you do, please consider 
>> providing a patch to the deal.II development sources (see the deal.II 
>> website on how to contribute).
>>
>> Stacktrace:
>> ---
>> #0  /mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre: 
>> dealii::MatrixTools::apply_boundary_values(std::map> std::complex, std::less, 
>> std::allocator > > > 
>> const&, dealii::PETScWrappers::MatrixBase&, 
>> dealii::PETScWrappers::VectorBase&, dealii::PETScWrappers::VectorBase&, 
>> bool)
>> #1  ./pFEM: Adaptive::LaplaceProblem<2>::assemble_system()
>> #2  ./pFEM: Adaptive::LaplaceProblem<2>::run()
>> #3  ./pFEM: main
>> 
>>
>> [jc-umit:02419] *** Process received signal ***
>> [jc-umit:02419] Signal: Aborted (6)
>> [jc-umit:02419] Signal code:  (-6)
>> [jc-umit:02419] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x36cb0) 
>> [0x7f80d19efcb0]
>> [jc-umit:02419] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37) 
>> [0x7f80d19efc37]
>> [jc-umit:02419] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x148) 
>> [0x7f80d19f3028]
>> [jc-umit:02419] [ 3] 
>> /mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(+0x8047862) 
>> [0x7f80da96e862]
>> [jc-umit:02419] [ 4] 
>> /mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseE+0x1a)
>>  
>> [0x7f80da96fe81]
>> [jc-umit:02419] [ 5] 
>> /mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_18StandardExceptions17ExcNotImplementedEEEvNS1_17ExceptionHandlingEPKciS7_S7_S7_T_+0x4e)
>>  
>> [0x7f80d5454e7e]
>> [jc-umit:02419] [ 6] 
>> /mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii11MatrixTools21apply_boundary_valuesERKSt3mapIjSt7complexIdESt4lessIjESaISt4pairIKjS3_EEERNS_13PETScWrappers10MatrixBaseERNSD_10VectorBaseESH_b+0xc2)
>>  
>> [0x7f80d54cedf6]
>> [jc-umit:02419] [ 7] 
>> ./pFEM(_ZN8Adaptive14LaplaceProblemILi2EE15assemble_systemEv+0xb47) 
>> [0x43a137]
>> [jc-umit:02419] [ 8] 
>> ./pFEM(_ZN8Adaptive14LaplaceProblemILi2EE3runEv+0x126) [0x43ce44]
>> [jc-umit:02419] [ 9] ./pFEM(main+0xdf) [0x427610]
>> [jc-umit:02419] [10] 
>> /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5) [0x7f80d19daf45]
>> [jc-umit:02419] [11] ./pFEM

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Thanks again, 

Now it passes the *interpolate_boundary_values *call. However, it is 
conflicting now with:

  MatrixTools::apply_boundary_values (bval,

system_matrix,

solution,

system_rhs);

To be precise I write the output below. 

Thanks in advance,
Juan Carlos Araújo, Umeå University



An error occurred in line <70> of file 
 in function
void dealii::MatrixTools::apply_boundary_values(const std::map >&, dealii::PETScWrappers::MatrixBase&, 
dealii::PETScWrappers::VectorBase&, dealii::PETScWrappers::VectorBase&, 
bool)
The violated condition was: 
eliminate_columns == false
Additional information: 
You are trying to use functionality in deal.II that is currently not 
implemented. In many cases, this indicates that there simply didn't appear 
much of a need for it, or that the author of the original code did not have 
the time to implement a particular case. If you hit this exception, it is 
therefore worth the time to look into the code to find out whether you may 
be able to implement the missing functionality. If you do, please consider 
providing a patch to the deal.II development sources (see the deal.II 
website on how to contribute).

Stacktrace:
---
#0  /mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre: 
dealii::MatrixTools::apply_boundary_values(std::map, std::less, 
std::allocator > > > 
const&, dealii::PETScWrappers::MatrixBase&, 
dealii::PETScWrappers::VectorBase&, dealii::PETScWrappers::VectorBase&, 
bool)
#1  ./pFEM: Adaptive::LaplaceProblem<2>::assemble_system()
#2  ./pFEM: Adaptive::LaplaceProblem<2>::run()
#3  ./pFEM: main


[jc-umit:02419] *** Process received signal ***
[jc-umit:02419] Signal: Aborted (6)
[jc-umit:02419] Signal code:  (-6)
[jc-umit:02419] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x36cb0) 
[0x7f80d19efcb0]
[jc-umit:02419] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37) 
[0x7f80d19efc37]
[jc-umit:02419] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x148) 
[0x7f80d19f3028]
[jc-umit:02419] [ 3] 
/mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(+0x8047862) 
[0x7f80da96e862]
[jc-umit:02419] [ 4] 
/mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseE+0x1a)
 
[0x7f80da96fe81]
[jc-umit:02419] [ 5] 
/mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_18StandardExceptions17ExcNotImplementedEEEvNS1_17ExceptionHandlingEPKciS7_S7_S7_T_+0x4e)
 
[0x7f80d5454e7e]
[jc-umit:02419] [ 6] 
/mypath/Soft/dealii/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii11MatrixTools21apply_boundary_valuesERKSt3mapIjSt7complexIdESt4lessIjESaISt4pairIKjS3_EEERNS_13PETScWrappers10MatrixBaseERNSD_10VectorBaseESH_b+0xc2)
 
[0x7f80d54cedf6]
[jc-umit:02419] [ 7] 
./pFEM(_ZN8Adaptive14LaplaceProblemILi2EE15assemble_systemEv+0xb47) 
[0x43a137]
[jc-umit:02419] [ 8] ./pFEM(_ZN8Adaptive14LaplaceProblemILi2EE3runEv+0x126) 
[0x43ce44]
[jc-umit:02419] [ 9] ./pFEM(main+0xdf) [0x427610]
[jc-umit:02419] [10] 
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5) [0x7f80d19daf45]
[jc-umit:02419] [11] ./pFEM() [0x427419]
[jc-umit:02419] *** End of error message ***
Aborted (core dumped)




El viernes, 8 de diciembre de 2017, 13:07:28 (UTC-5), Wolfgang Bangerth 
escribió:
>
>
> > Merge: fd97710 cceeb7c 
> > Author: Wolfgang Bangerth  
> > Date:   Tue Dec 5 07:24:15 2017 -0700 
> > 
> >  Merge pull request #5579 from bangerth/more-complex-instantiations 
> > 
> >  More complex instantiations 
> > 
> > 
> > However, after configuration, installation and compilation of my code I 
> > obtain the error: 
> > 
> > Scanning dependencies of target pFEM 
> > [100%] Building CXX object CMakeFiles/pFEM.dir/pFEM.cc.o 
> > Linking CXX executable pFEM 
> > /mypath/pFEM.cc:370: error: undefined reference to 'void 
> > dealii::VectorTools::interpolate_boundary_values<2, 2, 
> > std::complex >(dealii::hp::MappingCollection<2, 2> const&, 
> > dealii::hp::DoFHandler<2, 2> const&, std::map > dealii::Function<2, std::complex > const*, std::less > int>, std::allocator > std::complex > const*> > > const&, std::map > std::complex, std::less, 
> > std::allocator > > 
> >  >&, dealii::ComponentMask const&)' 
>
> Yes, I apparently forgot to instantiate this one function as well. Patch 
> is here: 
>https://github.com/dealii/dealii/pull/5602 
> It was just merged. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Hi, when I do
git log
the #5579 does show up. 

Merge: fd97710 cceeb7c
Author: Wolfgang Bangerth 
Date:   Tue Dec 5 07:24:15 2017 -0700

Merge pull request #5579 from bangerth/more-complex-instantiations

More complex instantiations


However, after configuration, installation and compilation of my code I 
obtain the error:

Scanning dependencies of target pFEM
[100%] Building CXX object CMakeFiles/pFEM.dir/pFEM.cc.o
Linking CXX executable pFEM
/mypath/pFEM.cc:370: error: undefined reference to 'void 
dealii::VectorTools::interpolate_boundary_values<2, 2, std::complex 
>(dealii::hp::MappingCollection<2, 2> const&, dealii::hp::DoFHandler<2, 2> 
const&, std::map > 
const*, std::less, std::allocator > const*> > > const&, 
std::map, std::less, 
std::allocator > > >&, 
dealii::ComponentMask const&)'
collect2: error: ld returned 1 exit status
make[2]: *** [pFEM] Error 1
make[1]: *** [CMakeFiles/pFEM.dir/all] Error 2
make: *** [all] Error 2

Thanks in advance,
Juan Carlos Araújo


El martes, 5 de diciembre de 2017, 23:47:22 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 12/05/2017 05:47 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > On my topic, I have done as usual: 
> >  git clone https://github.com/dealii/dealii.git dealii 
> > and proceeded with the configuration/installation of the library, to 
> notice 
> > that I still have an error: 
>
> When did you clone the repository? The patch was merged earlier today. You 
> can 
> find out if you do 
>git log 
> and look at which patches were merged most recently. Does #5579 show up? 
>
> It's of course also possible that I fixed one issue but not all. 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Thanks for your email and the encouragement, I will keep these reporting 
issues for sure!
This thread was a nice lecture on templates!

On my topic, I have done as usual:
git clone https://github.com/dealii/dealii.git dealii
and proceeded with the configuration/installation of the library, to notice 
that I still have an error:

Linking CXX executable pFEM
/mypatch/pFEM.cc:370: error: undefined reference to 'void 
dealii::VectorTools::interpolate_boundary_values<2, 2, std::complex 
>(dealii::hp::MappingCollection<2, 2> const&, dealii::hp::DoFHandler<2, 2> 
const&, std::map > 
const*, std::less, std::allocator > const*> > > const&, 
std::map, std::less, 
std::allocator > > >&, 
dealii::ComponentMask const&)'
collect2: error: ld returned 1 exit status
make[2]: *** [pFEM] Error 1
make[1]: *** [CMakeFiles/pFEM.dir/all] Error 2
make: *** [all] Error 2

As if nothing had changed. So I have a few questions:
1) Is it the same error?
2) If 1) is true, then isn't https://github.com/dealii/dealii/pull/5579 
<https://www.google.com/url?q=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fpull%2F5579&sa=D&sntz=1&usg=AFQjCNFmtOPSaCSooQRONvwAcX3PeaAI4A>
 
merged with https://github.com/dealii/dealii.git dealii?
3) If 2) is false, I would like to kindly ask how to apply the patch. I 
tried a few things I got from a GIT tutorial but seemed not to work.

Thanks in advance,
Juan Carlos Araújo, Umeå Universitet


El lunes, 4 de diciembre de 2017, 13:09:28 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 12/01/2017 08:36 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > path/petsc_eigs/pFEM.cc:370: error: undefined reference to 'void 
> > dealii::VectorTools::interpolate_boundary_values<2, 2, 
> > std::complex >(dealii::hp::MappingCollection<2, 2> const&, 
> > dealii::hp::DoFHandler<2, 2> const&, std::map > dealii::Function<2, std::complex > const*, std::less > int>, std::allocator > std::complex > const*> > > const&, std::map > std::complex, std::less, 
> > std::allocator > > 
> >  >&, dealii::ComponentMask const&)' 
>
> Now that the compiler is happy with your code, this is the first real 
> bug -- we don't instantiate this function. The patch is here: 
>https://github.com/dealii/dealii/pull/5579 
>
> I appreciate you trying all of these things with complex-valued vectors. 
> It's hard to find all of the relevant places that need to be changed and 
> fixed. Keep these issues coming! 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Thanks again! that one was very logical ... silly me!

To review, now it looks like this:

  std::map bval;
  
  ZeroFunction   homogeneous_dirichlet_bc;
  
  const typename FunctionMap::type 
dirichlet_boundary_functions
  = { { types::boundary_id(0), &homogeneous_dirichlet_bc } };   
   
  
  VectorTools::interpolate_boundary_values (mapping_collection,

  dof_handler,

dirichlet_boundary_functions,
bval);

and gives the following error (no candidates this time):

path/petsc_eigs/pFEM.cc:370: error: undefined reference to 'void 
dealii::VectorTools::interpolate_boundary_values<2, 2, std::complex 
>(dealii::hp::MappingCollection<2, 2> const&, dealii::hp::DoFHandler<2, 2> 
const&, std::map > 
const*, std::less, std::allocator > const*> > > const&, 
std::map, std::less, 
std::allocator > > >&, 
dealii::ComponentMask const&)'
collect2: error: ld returned 1 exit status
make[2]: *** [pFEM] Error 1
make[1]: *** [CMakeFiles/pFEM.dir/all] Error 2
make: *** [all] Error 2

Any comments? thanks in advance,
Juan Araújo


El viernes, 1 de diciembre de 2017, 15:49:48 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 12/01/2017 01:29 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > I must admit that I am a bit lost here. However, I tried your suggestion 
> > and used:FunctionMap::type,  what gave the errors 
> > that I attach in errors.txt. 
> > 
> > By looking at those errors, I also tried: 
> >  std::map > std::complex > >  bval; 
> > but it seems that this is not the root of the problem (I also attach 
> > these errors at the end of the same file). 
> > 
> > Any ideas? Thanks in advance. 
>
> Yes. You are now trying assign a `ZeroFunction` object's address to 
> a FunctionMap::type, but this wants the address of a 
> complex-valued function. You need to use 
>ZeroFunction 
> now to explain to the compiler that you want a function that returns a 
> std::complex (==PetscScalar)-typed zero, rather than a real-valued one. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Apologies for the last email ...

Thanks W. for your quick reply.

I used static because it was used in one of the dealii steps. As there is 
no need for it, in the following it is removed.

I must admit that I am a bit lost here. However, I tried your suggestion 
and used:FunctionMap::type,  what gave the errors that 
I attach in errors.txt.

By looking at those errors, I also tried:
std::map > >  bval;
but it seems that this is not the root of the problem (I also attach these 
errors at the end of the same file).

Any ideas? Thanks in advance.

El viernes, 1 de diciembre de 2017, 15:06:29 (UTC-5), Bruno Turcksin 
escribió:
>
>
> On Friday, December 1, 2017 at 1:47:44 PM UTC-5, Juan Carlos Araujo 
> Cabarcas wrote:
>
>> I must admit that I am a bit lost here. However, I tried your suggestion 
>> and used:FunctionMap::type,  what gave the errors that 
>> I attach in errors.txt.
>>
>> You have attached a file but it only contains an e-mail address. 
>
> Best,
>
> Bruno
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
** Errors using:  FunctionMap::type

jc/adapt/petsc_eigs/pFEM.cc: In instantiation of ‘void 
Adaptive::LaplaceProblem::assemble_system() [with int dim = 2]’:
jc/adapt/petsc_eigs/pFEM.cc:562:26:   required from ‘void 
Adaptive::LaplaceProblem::run() [with int dim = 2]’
jc/adapt/petsc_eigs/pFEM.cc:595:28:   required from here
jc/adapt/petsc_eigs/pFEM.cc:368:72: error: could not convert ‘{{0u, (& 
homogeneous_dirichlet_bc)}}’ from ‘’ to ‘const 
type {aka const std::map >*, std::less, 
std::allocator >*> > >}’
   = { { types::boundary_id(0), &homogeneous_dirichlet_bc } };
^
In file included from dealii/include/deal.II/base/point.h:21:0,
 from dealii/include/deal.II/base/quadrature.h:21,
 from dealii/include/deal.II/base/quadrature_lib.h:21,
 from jc/adapt/petsc_eigs/pFEM.cc:29:
dealii/include/deal.II/lac/constraint_matrix.h: In instantiation of ‘void 
dealii::ConstraintMatrix::distribute_local_to_global(const InVector&, const 
std::vector&, OutVector&) const [with InVector = 
dealii::FullMatrix; OutVector = dealii::PETScWrappers::SparseMatrix]’:
jc/adapt/petsc_eigs/pFEM.cc:352:50:   required from ‘void 
Adaptive::LaplaceProblem::assemble_system() [with int dim = 2]’
jc/adapt/petsc_eigs/pFEM.cc:562:26:   required from ‘void 
Adaptive::LaplaceProblem::run() [with int dim = 2]’
jc/adapt/petsc_eigs/pFEM.cc:595:28:   required from here
dealii/include/deal.II/lac/constraint_matrix.h:1761:31: error: no match for 
‘operator==’ (operand types are ‘const dealii::TableIndices<2>’ and 
‘std::vector::size_type {aka long unsigned int}’)
   Assert (local_vector.size() == local_dof_indices.size(),
   ^
dealii/include/deal.II/base/exceptions.h:339:11: note: in definition of macro 
‘Assert’
 if (!(cond))  \
   ^
dealii/include/deal.II/lac/constraint_matrix.h:1761:31: note: candidates are:
   Assert (local_vector.size() == local_dof_indices.size(),
   ^
dealii/include/deal.II/base/exceptions.h:339:11: note: in definition of macro 
‘Assert’
 if (!(cond))  \
   ^
In file included from dealii/include/deal.II/base/tensor.h:22:0,
 from dealii/include/deal.II/base/point.h:22,
 from dealii/include/deal.II/base/quadrature.h:21,
 from dealii/include/deal.II/base/quadrature_lib.h:21,
 from jc/adapt/petsc_eigs/pFEM.cc:29:
dealii/include/deal.II/base/table_indices.h:370:1: note: bool 
dealii::TableIndices::operator==(const dealii::TableIndices&) const [with 
int N = 2]
 TableIndices::operator == (const TableIndices &other) const
 ^
dealii/include/deal.II/base/table_indices.h:370:1: note:   no known conversion 
for argument 1 from ‘std::vector::size_type {aka long unsigned 
int}’ to ‘const dealii::TableIndices<2>&’
In file included from dealii/include/deal.II/fe/mapping_q_generic.h:24:0,
 from dealii/include/deal.II/fe/mapping_q1.h:21,
 from dealii/include/deal.II/hp/mapping_collection.h:21,
 from dealii/include/deal.II/numerics/vector_tools.h:27,
 from jc/adapt/petsc_eigs/pFEM.cc:61:
dealii/include/deal.II/base/vectorization.h:2719:1: note: tem

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Thanks for your quick reply.

I used static because it was used in one of the dealii steps. As there is 
no need for it, in the following it is removed.

I must admit that I am a bit lost here. However, I tried your suggestion 
and used:FunctionMap::type,  what gave the errors that 
I attach in errors.txt.

By looking at those errors, I also tried:
std::map > >  bval;
but it seems that this is not the root of the problem (I also attach these 
errors at the end of the same file).

Any ideas? Thanks in advance.

El viernes, 1 de diciembre de 2017, 8:09:26 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 11/30/2017 03:05 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > When using the PETSc wrappers as in step-36 (installed with PetscScalar 
> = 
> > complex) and using MappingCollection as discussed before. I now 
> try: 
> > 
> > | 
> > //static std::map bval; 
> > staticstd::mapbval; 
>
> This is unrelated to the error, but why do you make this variable 
> `static`? 
>
> > ZeroFunction   homogeneous_dirichlet_bc; 
> > 
> > consttypenameFunctionMap::type dirichlet_boundary_functions 
> > ={{types::boundary_id(0),&homogeneous_dirichlet_bc }}; 
>
> If you use a map for `bval` that has complex-valued keys, then you also 
> need 
> to use a FunctionMap::type object here so that the 
> boundary 
> functions are complex-valued as well. 
>
> In other words, the scalar type of the function object and of the boundary 
> values must match. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
andres.alf...@correounivalle.edu.co

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Thanks for your quick reply.

I used static because it was used like that in one of the steps, as there 
is no need for it, in the following it is removed.

I must admit that I am a bit lost here. However, I tried your suggestion 
and used:FunctionMap::type,  what gave the errors that 
I attach in errors.txt.
By looking at those errors, I also tried:

  std::map > >  bval;

but it seems that is not the root of the problem (I also attach these 
errors at the end of the same file).

Any ideas? Thanks in advance.


El viernes, 1 de diciembre de 2017, 8:09:26 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 11/30/2017 03:05 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > When using the PETSc wrappers as in step-36 (installed with PetscScalar 
> = 
> > complex) and using MappingCollection as discussed before. I now 
> try: 
> > 
> > | 
> > //static std::map bval; 
> > staticstd::mapbval; 
>
> This is unrelated to the error, but why do you make this variable 
> `static`? 
>
> > ZeroFunction   homogeneous_dirichlet_bc; 
> > 
> > consttypenameFunctionMap::type dirichlet_boundary_functions 
> > ={{types::boundary_id(0),&homogeneous_dirichlet_bc }}; 
>
> If you use a map for `bval` that has complex-valued keys, then you also 
> need 
> to use a FunctionMap::type object here so that the 
> boundary 
> functions are complex-valued as well. 
>
> In other words, the scalar type of the function object and of the boundary 
> values must match. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

May the same bug exist for the complex arithmetic case?

When using the PETSc wrappers as in step-36 (installed with PetscScalar = 
complex) and using MappingCollection as discussed before. I now try:

  //static std::map bval;
  static std::map bval;
  ZeroFunctionhomogeneous_dirichlet_bc;
  
  const typename FunctionMap::type dirichlet_boundary_functions
  = { { types::boundary_id(0), &homogeneous_dirichlet_bc } };   
   
  
  VectorTools::interpolate_boundary_values (mapping_collection,

  dof_handler,

dirichlet_boundary_functions,//typename FunctionMap::type(),
bval);
and obtain the errors I leave at the end. The most relevant line is:

numerics/vector_tools.h:1006:3: note:   template argument 
deduction/substitution failed:
pFEM.cc:374:51: note:   deduced conflicting types for parameter ‘number’ 
(‘double’ and ‘std::complex’)
   bval);
   ^
Errors: in attached file.

Any ideas?
Thanks in advance.







El domingo, 12 de noviembre de 2017, 18:32:00 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 11/10/2017 08:15 AM, Juan Carlos Araujo Cabarcas wrote: 
> > Yes, it explains the Function map better, and also thanks for the 
> pointer to 
> > step-16. Unfortunately, I still cannot get it to run =( ... By looking 
> at 
> > step-35, I got the line: 
> > | 
> > staticstd::mapbval; 
> > | 
> > 
> > With this, I now write: 
> > | 
> > ZeroFunctionhomogeneous_dirichlet_bc; 
> > 
> > consttypenameFunctionMap::type dirichlet_boundary_functions 
> > ={{types::boundary_id(0),&homogeneous_dirichlet_bc }}; 
> > 
> > staticstd::mapbval; 
> > 
> > VectorTools::interpolate_boundary_values ( 
> >mapping_collection, 
> >dof_handler, 
> >dirichlet_boundary_functions,//typename 
> FunctionMap::type(), 
> >bval 
> > ); 
> > 
> > | 
> > 
> >  From which I get the error (without any candidates): 
> > 
> > error: undefined reference to 
> > 
> > 'void dealii::VectorTools::interpolate_boundary_values<2, 2, double>( 
> >  dealii::hp::MappingCollection<2, 2> const&, 
> >  dealii::hp::DoFHandler<2, 2> const&, 
> >  std::map< 
> >  unsigned char, 
> >  dealii::Function< 
> >  2, double> const*, 
> >  std::less, 
> >  std::allocator > dealii::Function<2, double> const*> > 
> >  > const&, 
> >  std::map< 
> >  unsigned int, 
> >  double, 
> >  std::less, 
> >  std::allocator > >&, 
> >  dealii::ComponentMask const& 
> > )' 
>
> Yes, that's a bug. The fix is here: 
>https://github.com/dealii/dealii/pull/5448 
> Thanks for reporting this! 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/home/ju4nk4/jc/codes/adaptivity/petsc_eigs/pFEM.cc: In instantiation of ‘void 
Adaptive::LaplaceProblem::assemble_system() [with int dim = 2]’:
/home/ju4nk4/jc/codes/adaptivity/petsc_eigs/pFEM.cc:559:26:   required from 
‘void Adaptive::LaplaceProblem::run() [with int dim = 2]’
/home/ju4nk4/jc/codes/adaptivity/petsc_eigs/pFEM.cc:592:28:   required from here
/home/ju4nk4/jc/codes/adaptivity/petsc_eigs/pFEM.cc:370:51: error: no matching 
function for call to 
‘interpolate_boundary_values(dealii::hp::MappingCollection<2, 2>&, 
dealii::hp::DoFHandler<2, 2>&, const type&, std::map >&)’
   bval);
   ^
/home/ju4nk4/jc/codes/adaptivity/petsc_eigs/pFEM.cc:370:51: note: cand

Re: [deal.II] Using DataOut with MappingCollection

[Juan Carlos Araujo Cabarcas]

Thanks for your quick reply and all the support on this thread. Programming 
with templates it's not really my cup of tea!

El martes, 28 de noviembre de 2017, 0:10:15 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 11/27/2017 08:13 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > /home/ju4nk4/Soft/dealii/include/deal.II/numerics/data_out.h:285:16: 
> note: 
> > void dealii::DataOut::build_patches(const 
> > dealii::Mapping > space_dimension>&, unsigned int, dealii::DataOut > DoFHandlerType>::CurvedCellRegion) [with int dim = 2; DoFHandlerType = 
> > dealii::hp::DoFHandler<2, 2>] 
> > virtual void build_patches (const Mapping > DoFHandlerType::space_dimension> &mapping, 
> >  ^ 
> > /home/ju4nk4/Soft/dealii/include/deal.II/numerics/data_out.h:285:16: 
> note:   
> > no known conversion for argument 3 from ‘dealii::DataOut<2, 
> > dealii::DoFHandler<2, 2> >::CurvedCellRegion’ to ‘dealii::DataOut<2, 
> > dealii::hp::DoFHandler<2, 2> >::CurvedCellRegion’ 
>
> Here's your key -- you need to say 
>DataOut >::curved_inner_cells 
> instead of 
>DataOut::curved_inner_cells 
>
> Error messages are your friend :-) 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Using DataOut with MappingCollection

[Juan Carlos Araujo Cabarcas]

Thanks for your reply. I did not know I could just pass the higher order 
mapping. If I omit DataOut::curved_inner_cells, then it works. 
However, I pursue domains with curved inner cells. Following your 
suggestion, I do:

data_out.build_patches (MappingQGeneric(max_degree), 8, DataOut::
curved_inner_cells);

and get the errors:

/home/ju4nk4/jc/codes/adaptivity/disk_eigs/pFEM.cc: In instantiation of 
‘void Adaptive::LaplaceProblem::postprocess(unsigned int) [with int 
dim = 2]’:
/home/ju4nk4/jc/codes/adaptivity/disk_eigs/pFEM.cc:502:27:   required from 
‘void Adaptive::LaplaceProblem::run() [with int dim = 2]’
/home/ju4nk4/jc/codes/adaptivity/disk_eigs/pFEM.cc:529:28:   required from 
here
/home/ju4nk4/jc/codes/adaptivity/disk_eigs/pFEM.cc:429:7: error: no 
matching function for call to ‘dealii::DataOut<2, dealii::hp::DoFHandler<2, 
2> >::build_patches(dealii::MappingQGeneric<2, 2>, int, dealii::DataOut<2, 
dealii::DoFHandler<2, 2> >::CurvedCellRegion)’
   data_out.build_patches (MappingQGeneric(max_degree), 8, 
DataOut::curved_inner_cells);
   ^
/home/ju4nk4/jc/codes/adaptivity/disk_eigs/pFEM.cc:429:7: note: candidates 
are:
In file included from 
/home/ju4nk4/jc/codes/adaptivity/disk_eigs/pFEM.cc:58:0:
/home/ju4nk4/Soft/dealii/include/deal.II/numerics/data_out.h:252:16: note: 
void dealii::DataOut::build_patches(unsigned int) 
[with int dim = 2; DoFHandlerType = dealii::hp::DoFHandler<2, 2>]
   virtual void build_patches (const unsigned int n_subdivisions = 0);
^
/home/ju4nk4/Soft/dealii/include/deal.II/numerics/data_out.h:252:16: 
note:   candidate expects 1 argument, 3 provided
/home/ju4nk4/Soft/dealii/include/deal.II/numerics/data_out.h:285:16: note: 
void dealii::DataOut::build_patches(const 
dealii::Mapping&, unsigned int, dealii::DataOut::CurvedCellRegion) [with int dim = 2; DoFHandlerType = 
dealii::hp::DoFHandler<2, 2>]
   virtual void build_patches (const Mapping &mapping,
^
/home/ju4nk4/Soft/dealii/include/deal.II/numerics/data_out.h:285:16: 
note:   no known conversion for argument 3 from ‘dealii::DataOut<2, 
dealii::DoFHandler<2, 2> >::CurvedCellRegion’ to ‘dealii::DataOut<2, 
dealii::hp::DoFHandler<2, 2> >::CurvedCellRegion’
make[2]: *** [CMakeFiles/pFEM.dir/pFEM.cc.o] Error 1
make[1]: *** [CMakeFiles/pFEM.dir/all] Error 2
make: *** [all] Error 2

Any ideas?


El lunes, 27 de noviembre de 2017, 12:45:04 (UTC-5), Wolfgang Bangerth 
escribió:
>
> On 11/17/2017 12:37 PM, Juan Carlos Araujo Cabarcas wrote: 
> > 
> > I would like to reproduce step-27 but with curved boundaries with the 
> > use of MappingCollection. 
> > Everything seems to work fine, but I noticed that data_out does not seem 
> > to be implemented for passing a MappingCollection. 
>
> Yes, that seems to be correct. There is even a @todo in the 
> documentation of the function that takes a mapping. 
>
>
> > In particular I would like to be able to use something like: 
> >data_out.build_patches (mapping_collection, 8, 
> > DataOut::curved_inner_cells); 
> > 
> > Any hints on how to achieve this are greatly appreciated! 
>
> I suspect that -- unless you are on a very coarse mesh -- the difference 
> between the different mappings are not really visible in a 
> visualization. Could you just pass the higher order mapping to the 
> function, instead of the entire mapping collection? 
>
> That's not the "correct" approach, of course, but it's likely not going 
> to lead to visible differences. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Using DataOut with MappingCollection

[Juan Carlos Araujo Cabarcas]

Dear all, 

As described in the forum:
https://groups.google.com/forum/?utm_source=digest&utm_medium=email#!topic/dealii/ijJFbOTZxkk

I would like to reproduce step-27 but with curved boundaries with the use 
of MappingCollection.
Everything seems to work fine, but I noticed that data_out does not seem to 
be implemented for passing a MappingCollection.

In particular I would like to be able to use something like:
  data_out.build_patches (mapping_collection, 8, 
DataOut::curved_inner_cells);

Any hints on how to achieve this are greatly appreciated!

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Thanks for the quick answer!
My problem is the scalar wave equation and the boundary conditions are 
simple homogeneous Dirichlet. Before starting with hp I was successfully 
using ZeroFunction(). How would this "*map taking a types::boundary_id 
and a Function object*" look like? 
Any step in the tutorial showing this?

Thanks,

El lunes, 23 de octubre de 2017, 15:08:27 (UTC-5), Daniel Arndt escribió:
>
> Juan
>
> [...]
>>  
>>static std::map bval;
>> VectorTools::interpolate_boundary_values (mapping_collection, 
>> dof_handler, ZeroFunction(), bval);
>>
>> and I get the following error:
>>
>> [100%] Building CXX object CMakeFiles/pFEM.dir/pFEM.cc.o
>> /home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc: In instantiation of 
>> ‘void Adaptive::LaplaceProblem::setup_system() [with int dim = 2]’:
>> /home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:377:23:   required 
>> from ‘void Adaptive::LaplaceProblem::run() [with int dim = 2]’
>> /home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:418:28:   required 
>> from here
>> /home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:205:51: error: no 
>> matching function for call to 
>> ‘interpolate_boundary_values(dealii::hp::MappingCollection<2, 2>&, 
>> dealii::hp::DoFHandler<2, 2>&, dealii::ZeroFunction<2, double>, 
>> std::map&)’
>>bval);
>>^
>> /home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:205:51: note: 
>> candidates are: ...
>>
> The compiler is your friend here. If you look at the possible candidates, 
> it tells why each of them could not be used.
>
> In particular, the signature for the function you are looking for is
>
> VectorTools::interpolate_boundary_values 
> (const hp::MappingCollection< dim, spacedim > &   
>   mapping,
>  const hp::DoFHandler< dim, spacedim > & 
> dof,
>  const std::map< types::boundary_id, const Function< spacedim, number > * 
> > & function_map,
>  std::map< types::global_dof_index, number > &   
>   boundary_values,
>  const ComponentMask &   
>   component_mask = ComponentMask())
>
> and the third parameter is a not a Function, but a map taking a 
> types::boundary_id and a Function object.
>
> Best,
> Daniel
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] mapping_collection and step-27

[Juan Carlos Araujo Cabarcas]

Dear all,

I would like to reproduce step-27 for a curved domain instead of the given 
domain.
For this, I would like to make use of a MappingCollection. 

I am having problems using interpolate_boundary_values, suggested in:
https://www.dealii.org/8.5.0/doxygen/deal.II/namespaceVectorTools.html#a074cd2642b4abdfd67db2c672c907781

Basically I use the following:

hp::DoFHandler  dof_handler;
hp::FECollectionfe_collection;
hp::MappingCollectionmapping_collection;
hp::QCollection quadrature_collection;
hp::QCollection   face_quadrature_collection;



and initialize them as step-27. Then I try:

 
   static std::map bval;
VectorTools::interpolate_boundary_values (mapping_collection, 
dof_handler, ZeroFunction(), bval);



and I get the following error:

[100%] Building CXX object CMakeFiles/pFEM.dir/pFEM.cc.o
/home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc: In instantiation of 
‘void Adaptive::LaplaceProblem::setup_system() [with int dim = 2]’:
/home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:377:23:   required from 
‘void Adaptive::LaplaceProblem::run() [with int dim = 2]’
/home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:418:28:   required from 
here
/home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:205:51: error: no 
matching function for call to 
‘interpolate_boundary_values(dealii::hp::MappingCollection<2, 2>&, 
dealii::hp::DoFHandler<2, 2>&, dealii::ZeroFunction<2, double>, 
std::map&)’
   bval);
   ^
/home/ju4nk4/jc/codes/adaptivity/disk_pFEM/pFEM.cc:205:51: note: candidates 
are: ...

I digged into the documentation for a while before emailing, but found that 
MappingCollection is not really discussed much, and the documentation for 
the call VectorTools::interpolate_boundary_values for MappingCollection 
just says:

"Like the previous function, but take a mapping collection to go with 
the hp::DoFHandler object." ...

So, there is no special comments there either.

Any guidance would be very much appreciated!

Juan Carlos Araújo,
Umeå Universitet

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] lexicographic ordering in 3D

[Juan Carlos Araujo Cabarcas]

Thanks for your quick and useful answer, and explanation ... GridReordering 
was what I needed!
Now I see that I did not mean lexicographic order, that I had assumed from 
the begining in 2), but it didn't read like that in the email.
I apologize for the wrong formulation of the question.

*Some constructive criticism though* ... In the documentation of 
GridReordering,
https://www.dealii.org/8.4.0/doxygen/deal.II/classGridReordering.html

The statement: ... "Hence, for creating a Triangulation 
 
based on the resulting CellData 
 the 
Triangulation::create_triangulation_compatibility() 

 
(and not the Triangulation::create_triangulation() 
)
 
function must be used."

should be stated more times and in places where it would be easier to 
reach, like in the description of GridReordering::reorder_cells() 

 
in the *Member Function Documentation*, or/and in the exceptions of 
Triangulation::create_triangulation() 
.
  
It took me a while to realize that it had to be done!

El lunes, 10 de julio de 2017, 19:45:45 (UTC+2), Wolfgang Bangerth escribió:
>
>
> Juan Carlos, 
>
> > The procedure goes as follows, 
> > 1) we introduce a set of vertices as points (x_i,y_i), numbered by the 
> > index "i", 
> > 2) connect them forming quads/cells and store the numbering as: 
> > cell_nodes [k] = {i1,i2,i3,i4}, 
> > 3) making sure that i1,i2,i3,i4 follow a lexicographic ordering 
> > explained in .../structGeometryInfo.html. 
> > 
> > In 2D it is messy already to satisfy 1), 2), 3) for 
> > non-trivial/primitive meshes, but we always can sketch the situation, 
> > because it is 2D, and by rotating cells (ordering) we hope that 
> > eventually a solution will come out. 
> > 
> > Now, the sketch part is not available in 3D and the number of 
> > permutations/rotations has increased dramatically. This leads me to the 
> > following question: 
> > 
> > Is there a tool within deal.II to impose 3), when 1) and 2) are given? 
> > In other words, is there any tool available, such that we get an array 
> > new_cell_nodes satisfying a lexicographic ordering, by inputting the 
> > arrays: vertex, cell_nodes? 
>
> In 2d, it would in principle be possible to write such code that, given 
> an unordered set of vertices brings them into an order so that they form 
> a lexicographic ordering. That's because -- at least for convex 
> quadrilaterals -- you can find the centroid of the vertices, sort them 
> in counter-clockwise sense, and then swap the last two. 
>
> But that can't work in 3d. Think, for example, of the vertices of a 
> cube. Number them as you usually would, and you'd get a cube -- a valid 
> geometry. But then rotate the top surface by 90 degrees: the (unsorted 
> set of) vertices are in the exact same location, but the geometry you 
> want to describe (a "twisted" cube) is a different, though equally 
> valid, one. In other words, you can't infer the geometry someone has in 
> mind just from knowing where the vertices are. Consequently, it is not 
> possible to write a code that can always get that right. 
>
>
> The only thing we can do for you is to avoid having to deal with the 
> tedium of orienting cells relative to each other. This is what the 
> GridReordering class does for you. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] lexicographic ordering in 3D

[Juan Carlos Araujo Cabarcas]

Dear all,

I have been designing 2D grids according to the standards given in:
https://www.dealii.org/8.4.1/doxygen/deal.II/structGeometryInfo.html

The procedure goes as follows,
1) we introduce a set of vertices as points (x_i,y_i), numbered by the 
index "i",
2) connect them forming quads/cells and store the numbering as: cell_nodes 
[k] = {i1,i2,i3,i4},
3) making sure that i1,i2,i3,i4 follow a lexicographic ordering explained 
in .../structGeometryInfo.html.

In 2D it is messy already to satisfy 1), 2), 3) for non-trivial/primitive 
meshes, but we always can sketch the situation, because it is 2D, and by 
rotating cells (ordering) we hope that eventually a solution will come out.

Now, the sketch part is not available in 3D and the number of 
permutations/rotations has increased dramatically. This leads me to the 
following question:

Is there a tool within deal.II to impose 3), when 1) and 2) are given? 
In other words, is there any tool available, such that we get an array 
new_cell_nodes satisfying a lexicographic ordering, by inputting the 
arrays: vertex, cell_nodes?

Ideas on how to achieve this are greatly appreciated!

Thanks in advance,
Juan Carlos Araújo, Umeå Universitet

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to initialize PETScWrappers::SparseMatrix from a PETSc Mat

[Juan Carlos Araujo Cabarcas]

Thanks a lot Praveen, this is what I needed!

El martes, 6 de junio de 2017, 14:17:02 (UTC+2), Praveen C escribió:
>
> In my applications, I add this to CMakeLists.txt
>
> # In release mode, disable some warning messages
>
> IF(CMAKE_BUILD_TYPE MATCHES Release)
>
> SET( CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Wno-tautological-compare" )
>
> SET( CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Wno-unused-parameter" )
>
> SET( CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Wno-unused-variable" )
>
> ENDIF()
>
> Then I will get warnings in Debug mode but not in Release mode.
>
> praveen
>
> On Tue, Jun 6, 2017 at 8:07 AM, Bruno Turcksin  > wrote:
>
>> Juan Carlos
>>
>> On Tuesday, June 6, 2017 at 5:34:09 AM UTC-4, Juan Carlos Araujo Cabarcas 
>> wrote:
>>>
>>> The warnings come from my code, which is in developing stage! 
>>> In the installation I had before, 
>>> [-Wunused-variable],[-Wunused-parameter],[-Wunused-but-set-variable] did 
>>> not show up. It's a bummer that changing deal.II version will make such a 
>>> difference! As I am not proficient in cmake, I still think it would be 
>>> useful to know how to pass compiler flags by modifying CMakeLists.
>>>
>> cmake -DCMAKE_CXX_FLAGS="my_flags" should work.
>>
>> Best,
>>
>> Bruno 
>>
>> -- 
>> 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 .
>> For more options, visit https://groups.google.com/d/optout.
>>
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to initialize PETScWrappers::SparseMatrix from a PETSc Mat

[Juan Carlos Araujo Cabarcas]

Thanks for answering!

The warnings come from my code, which is in developing stage! 
In the installation I had before, 
[-Wunused-variable],[-Wunused-parameter],[-Wunused-but-set-variable] did 
not show up. It's a bummer that changing deal.II version will make such a 
difference! As I am not proficient in cmake, I still think it would be 
useful to know how to pass compiler flags by modifying CMakeLists.

El lunes, 5 de junio de 2017, 22:48:28 (UTC+2), Wolfgang Bangerth escribió:
>
>
> Juan Carlos, 
>
> > 1) First, I get way too many warnings like: 
> >   
> [-Wunused-variable],[-Wunused-parameter],[-Wunused-but-set-variable],etc 
> ... 
> > 
> > Then I would like to ask if there is a clean to disable warnings of this 
> type? 
> > maybe adding some lines to the CMakeLists? 
>
> Are these warnings coming out of deal.II, or out of your own code? If the 
> latter, you should of course fix them -- warnings warn you about issues 
> with 
> your code for very good reasons! 
>
>
> > 2) I realized that the type *PETScWrappers::Vector *has disappeared from 
> the 
> > library! 
> > As mentioned above, I was comfortably reading a native PETSc vector 
> *vec* from 
> > disk by doing *PETScWrappers::Vector (vec); *Now that it is gone,**I 
> don't 
> > really know what is the successor. 
>
> Just use PETScWrappers::MPI::Vector, but with MPI_COMM_SELF as the 
> communicator. This creates a local vector. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to initialize PETScWrappers::SparseMatrix from a PETSc Mat

[Juan Carlos Araujo Cabarcas]

Well, now I tried to update my code ... I was using:
git clone https://github.com/davydden/dealii
cd dealii
git checkout branch_petscscalar_complex

Now, because I would like to use a few recent additions from the 
development version, I do instead:
git clone https://github.com/dealii/dealii

1) First, I get way too many warnings like:
 
[-Wunused-variable],[-Wunused-parameter],[-Wunused-but-set-variable],etc ...

Then I would like to ask if there is a clean to disable warnings of this 
type? maybe adding some lines to the CMakeLists?

2) I realized that the type *PETScWrappers::Vector *has disappeared from 
the library! 
As mentioned above, I was comfortably reading a native PETSc vector *vec* 
from disk by doing *PETScWrappers::Vector (vec); *Now that it is gone, I 
don't really know what is the successor. 

By looking in step-17, *PETScWrappers::MPI::Vector *is passed to Vector for 
serial computations:
Vector 
 
localized_solution (solution);
I guess I could code a function to read native PETSc vectors into 
 
Vector
 
*, * but I 
still think that if the latter is the successor of *PETScWrappers::Vector, 
*then 
it is intuitive to have an equivalent constructor like existed before in 
*PETScWrappers::Vector.* 



-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: MappingQ, and minimal test for reproducing Benchmark

[Juan Carlos Araujo Cabarcas]

Thanks a lot Martin, 

I had been waiting for this type of solution in dealii for a long time!
I also read in another thread about using MappingQGeneric instead of 
MappingQ, so I updated my test code as below, and everything worked like a 
charm!


PolarManifold boundary;
TransfiniteInterpolationManifold inner_manifold;

  Triangulation triangulation;
  
Geom_parameters gp;
unsigned int curved_faces_label=10, not_curved = 1;

cell_rod ( triangulation, 1.0, R, gp, false );

triangulation.set_manifold ( curved_faces_label, 
boundary ); // labels defined with the grid
  
  inner_manifold.initialize(triangulation);
  triangulation.set_manifold (not_curved, inner_manifold);
  
  const MappingQGeneric mapping (degree);
const QGauss quadrature(degree);
  
  // const FE_Q dummy_fe (1);
const FE_Q dummy_fe (QGaussLobatto<1>(degree + 1));
  DoFHandler dof_handler (triangulation);

  FEValues fe_values (mapping, dummy_fe, quadrature,
   update_JxW_values);




Juan Carlos Araújo,
Umeå Universitet

El sábado, 3 de junio de 2017, 10:36:58 (UTC+2), Martin Kronbichler 
escribió:
>
> Dear Juan Carlos,
>
> We have now finalized the TransfiniteInterpolationManifold implementation 
> in deal.II, 
>
> http://dealii.org/developer/doxygen/deal.II/classTransfiniteInterpolationManifold.html
>
> If you are interested, please try it out and let us know. As its 
> implementation is completely new and I have only applied it to two 
> applications in my group, there might still be issues, so let us know.
>
> By the way, if you are working on a 2D disk I recommend using 
> PolarManifold rather than SphericalManifold. I have not yet looked in your 
> code, but I recommend to try the transfinite interpolation manifold first.
>
> Best,
> Martin
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to initialize PETScWrappers::SparseMatrix from a PETSc Mat

[Juan Carlos Araujo Cabarcas]

Yeah, that was exactly my point, I just thought that it was already 
implemented but I just did not find the correct method or constructor.
I don't think the code I provided above uses all high performance 
capabilities available from PETSc ... I am in a testing phase so it is good 
enough for me.

El viernes, 26 de mayo de 2017, 19:51:56 (UTC+2), David Wells escribió:
>
> I think that this is a use case that we should support: for example, 
> PETScWrappers::VectorBase can act as a wrapper around a user-provided Vec 
> so we should definitely consider providing similar functionality for 
> PETScWrappers::MatrixBase.
>
> On Wednesday, May 24, 2017 at 5:18:46 AM UTC-4, Juan Carlos Araujo 
> Cabarcas wrote:
>>
>> Dear all,
>>
>> I have a matrix written by using PETSc in a Binary file, and I am trying 
>> work with it in my dealii environment. For this, I read a system matrix 
>> from PETSc:
>>
>> Mat   M_read;
>>
>> PetscViewerBinaryOpen(PETSC_COMM_WORLD,file.c_str (),FILE_MODE_READ,&fd);
>> MatCreate(PETSC_COMM_WORLD, &M_read);
>> MatSetOptionsPrefix( M_read,"m_");
>> MatSetFromOptions( M_read );
>> MatLoad( local_M ,fd);
>> PetscViewerDestroy(&fd);
>> MatGetSize( M_read ,&mm,&nn);
>> MatGetInfo( M_read ,MAT_LOCAL,&matinfo);
>> printf("matinfo.nz_used %g\n", matinfo.nz_used);
>>
>> and I would like to initialize a *dealii PETScWrappers::SparseMatrix *by 
>> using *M_read*.
>>
>> If I try the naĩve way:
>>   PETScWrappers::SparseMatrix M = PETScWrappers::SparseMatrix (M_read);
>>
>> I get the error:
>>
>> error: ‘dealii::PETScWrappers::SparseMatrix::SparseMatrix(const 
>> dealii::PETScWrappers::SparseMatrix&)’ is private
>>  SparseMatrix(const SparseMatrix &);
>>
>> I know the *PETScWrappers::SparseMatrix* contains the protected member *Mat 
>> matrix*, but I quite don't know how to copy it to initialize *M*.
>>
>> Is there a way to achieve this? any help is appreciated!
>>
>> Thanks in advance,
>> Juan Carlos Araújo,
>> Umeå Universitet
>>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to initialize PETScWrappers::SparseMatrix from a PETSc Mat

[Juan Carlos Araujo Cabarcas]

I managed to do it by looking at: *PETScWrappers::MatrixBase*, but my 
solution is extremely ugly!

  double max_coup_dofs = dof_handler.max_couplings_between_dofs();
  PETScWrappers::SparseMatrix local_M;
  local_M.reinit (par.dofs,par.dofs, max_coup_dofs);
  
  MatSetOption( local_M, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
  
{ // reading and parsing from (Mat M_read) into 
PETScWrappers::SparseMatrix local_M;
PetscInt begin, end;
const int ierr = MatGetOwnershipRange (
static_cast(M_read), &begin, &end);

std::pair
loc_range = std::make_pair (begin, end);

PetscInt ncols;
const PetscInt*colnums;
const PetscScalar *values;

PETScWrappers::MatrixBase::size_type row;
for (row = loc_range.first; row < loc_range.
second; ++row) {
 int ierr = MatGetRow(M_read, row, &ncols, &
colnums, &values);
 (void)ierr;
 AssertThrow (ierr == 0, PETScWrappers::
MatrixBase::ExcPETScError(ierr));
 
 for (PetscInt col = 0; col < ncols; ++col) 
{
 // cout << "(" << row << "," << 
colnums[col] << ") " << values[col] << std::endl;
 local_M.set (row,colnums[col], 
values[col] );
 }

 ierr = MatRestoreRow(M_read, row, &ncols, &
colnums, &values);
 AssertThrow (ierr == 0, PETScWrappers::
MatrixBase::ExcPETScError(ierr));
 }
//AssertThrow (out, ExcIO());
}

  local_M.compress (VectorOperation::add);
  //printf("\n\n\n  local_M norm is %f *** \n\n\n", 
local_M.l1_norm () );
  
  if (abs (result-local_M.l1_norm ()) > 1e-10) {
  printf("\n\n  M norm difference greater than 1e-10 after 
reading  *** \n\n" );
  } else {
  printf("\n\n Mass matrix successfully read! \n\n" );
  }

Any simpler alternative?

By the way, I noticed also the definition:
MatrixBase::operator Mat () const 
  { 
 return matrix; 
   }
By I was unable to use it!

Thanks,
Juan Carlos Araújo,
Umeå Universitet

El miércoles, 24 de mayo de 2017, 11:18:46 (UTC+2), Juan Carlos Araujo 
Cabarcas escribió:
>
> Dear all,
>
> I have a matrix written by using PETSc in a Binary file, and I am trying 
> work with it in my dealii environment. For this, I read a system matrix 
> from PETSc:
>
> Mat   M_read;
>
> PetscViewerBinaryOpen(PETSC_COMM_WORLD,file.c_str (),FILE_MODE_READ,&fd);
> MatCreate(PETSC_COMM_WORLD, &M_read);
> MatSetOptionsPrefix( M_read,"m_");
> MatSetFromOptions( M_read );
> MatLoad( local_M ,fd);
> PetscViewerDestroy(&fd);
> MatGetSize( M_read ,&mm,&nn);
> MatGetInfo( M_read ,MAT_LOCAL,&matinfo);
> printf("matinfo.nz_used %g\n", matinfo.nz_used);
>
> and I would like to initialize a *dealii PETScWrappers::SparseMatrix *by 
> using *M_read*.
>
> If I try the naĩve way:
>   PETScWrappers::SparseMatrix M = PETScWrappers::SparseMatrix (M_read);
>
> I get the error:
>
> error: ‘dealii::PETScWrappers::SparseMatrix::SparseMatrix(const 
> dealii::PETScWrappers::SparseMatrix&)’ is private
>  SparseMatrix(const SparseMatrix &);
>
> I know the *PETScWrappers::SparseMatrix* contains the protected member *Mat 
> matrix*, but I quite don't know how to copy it to initialize *M*.
>
> Is there a way to achieve this? any help is appreciated!
>
> Thanks in advance,
> Juan Carlos Araújo,
> Umeå Universitet
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] How to initialize PETScWrappers::SparseMatrix from a PETSc Mat

[Juan Carlos Araujo Cabarcas]

Dear all,

I have a matrix written by using PETSc in a Binary file, and I am trying 
work with it in my dealii environment. For this, I read a system matrix 
from PETSc:

Mat   M_read;

PetscViewerBinaryOpen(PETSC_COMM_WORLD,file.c_str (),FILE_MODE_READ,&fd);
MatCreate(PETSC_COMM_WORLD, &M_read);
MatSetOptionsPrefix( M_read,"m_");
MatSetFromOptions( M_read );
MatLoad( local_M ,fd);
PetscViewerDestroy(&fd);
MatGetSize( M_read ,&mm,&nn);
MatGetInfo( M_read ,MAT_LOCAL,&matinfo);
printf("matinfo.nz_used %g\n", matinfo.nz_used);

and I would like to initialize a *dealii PETScWrappers::SparseMatrix *by 
using *M_read*.

If I try the naĩve way:
  PETScWrappers::SparseMatrix M = PETScWrappers::SparseMatrix (M_read);

I get the error:

error: ‘dealii::PETScWrappers::SparseMatrix::SparseMatrix(const 
dealii::PETScWrappers::SparseMatrix&)’ is private
 SparseMatrix(const SparseMatrix &);

I know the *PETScWrappers::SparseMatrix* contains the protected member *Mat 
matrix*, but I quite don't know how to copy it to initialize *M*.

Is there a way to achieve this? any help is appreciated!

Thanks in advance,
Juan Carlos Araújo,
Umeå Universitet

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: state of dealii in complex arithmetics

[Juan Carlos Araujo Cabarcas]

Thanks, Denis and Wolfgang for the quick reply. 
I will give it a try and report back whatever issue I encounter.

El miércoles, 24 de mayo de 2017, 7:26:42 (UTC+2), Denis Davydov escribió:
>
>
> On Tuesday, May 23, 2017 at 5:22:44 PM UTC+2, Juan Carlos Araujo Cabarcas 
> wrote:
>>
>> Dear all,
>>
>> There are features in SLEPc 3.7 that I am interested in, and will try to 
>> make a fresh re-installation soon.
>>
>> I have been using deal.II 8.3 with complex arithmetic’s from the branch: 
>> git checkout branch_petscscalar_complex
>> and I wonder about the state of that branch, or if it has been merged to 
>> the master branch.
>>
>
> yes, everything in that branch is in the main repo from 8.4.2 or even 
> earlier.
>  
>
>>
>> Thanks in advance,
>>
>> Juan Carlos Araújo,
>> Umeå Universitet
>>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] state of dealii in complex arithmetics

[Juan Carlos Araujo Cabarcas]

Dear all,

There are features in SLEPc 3.7 that I am interested in, and will try to 
make a fresh re-installation soon.

I have been using deal.II 8.3 with complex arithmetic’s from the branch: 
git checkout branch_petscscalar_complex
and I wonder about the state of that branch, or if it has been merged to 
the master branch.

Thanks in advance,

Juan Carlos Araújo,
Umeå Universitet

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Arbitrary mapping or Manifold?

[Juan Carlos Araujo Cabarcas]

I wanted to add that ChartManifold along with MappingQ seems to be a very 
general way of parametrizing geometries and work for more general shapes 
than circles and ellipses (see attached .png). However, the corner points 
of hyper_ball proved to be very important at the time of performing the 
transformation. The figure I present has not been optimized (initial 
vertexes), but serves as motivation for illustrating the point.

I believe this thread is answered now, however I cannot mark it as answered 
... maybe an administrator can?

Thanks again,

El viernes, 7 de abril de 2017, 15:24:51 (UTC+2), Juan Carlos Araujo 
Cabarcas escribió:
>
> In the same spirit I would like to share a code that uses ChartManifold, 
> for generating a mapping of an ellipse applied to HyperBall.
>
> El miércoles, 5 de abril de 2017, 14:47:17 (UTC+2), Juan Carlos Araujo 
> Cabarcas escribió:
>>
>> Dear all, 
>> I recently found the thread: Something wrong with ChartManifold when 
>> chartdim=1:
>> https://groups.google.com/forum/#!topic/dealii/Sfo9xKoeRpw
>> a more straight answer to this thread.
>>
>> I took the liberty to copy David's code and modify it, so that the upper 
>> face of a square follows: y(x)=0.25*sin(PI*x) + 1. 
>> The function is not the same as I first proposed, but from this example 
>> it is pretty clear how to proceed. 
>> The result is plotted in the attached .png, along with a variation from 
>> step-10 in the deal.ii tutorial.
>>
>> I hope this may be of any use to somebody!
>>
>> Juan Carlos Araújo,
>> Umeå Universitet
>> El martes, 9 de febrero de 2016, 17:10:51 (UTC+1), Juan Carlos Araujo 
>> Cabarcas escribió:
>>>
>>> Dear all,
>>>
>>> I am working with the wave equation with variable refractive indexes 
>>> within the domain.
>>> I receive meshes that come from a shape optimization routine, and my 
>>> task is to
>>> run on arbitrarily curved elements resulting from the optimization.
>>> Step-10 and others show how to use mappings and Manifolds to curve 
>>> elements following basic shapes like circles, cylinders etc.
>>>
>>> My guess is that it should be possible in deal.II to define my own 
>>> Manifolds based on how I want to bend my edges (optimization routine).
>>>
>>> I have prepared a very basic example on the matlab code attached that 
>>> has the info to plot the attached figure. There, the edges follow the 
>>> equation |x|^p + |y|^p=1, with p=6.
>>>
>>> It would be very illustrative to apply a mapping like what is done in 
>>> step-10 on a circle, but for the presented case. I would appreciate any 
>>> advise on how to achieve this, hopefully wth this example =)
>>>
>>> Thanks in advance!
>>>
>>> Juan Carlos Araújo-Cabarcas
>>> Umeå Universitet.
>>>
>>>
>>>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Arbitrary mapping or Manifold?

[Juan Carlos Araujo Cabarcas]

I will see places where it may fit, I will come back to you.

Hej, Toby! nice to hear that it may help, looking forward to see the result 
=)



El viernes, 7 de abril de 2017, 15:42:13 (UTC+2), Wolfgang Bangerth 
escribió:
>
> On 04/07/2017 07:24 AM, Juan Carlos Araujo Cabarcas wrote: 
> > In the same spirit I would like to share a code that uses ChartManifold, 
> > for generating a mapping of an ellipse applied to HyperBall. 
>
> For both of these cases, do you think that there is a good place 
> somewhere in the documentation where this would fit? For the ellipse, 
> maybe as an example in the documentation of the ChartManifold class? 
> We'd be happy to add it there or somewhere else! 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Arbitrary mapping or Manifold?

[Juan Carlos Araujo Cabarcas]

In the same spirit I would like to share a code that uses ChartManifold, 
for generating a mapping of an ellipse applied to HyperBall.

El miércoles, 5 de abril de 2017, 14:47:17 (UTC+2), Juan Carlos Araujo 
Cabarcas escribió:
>
> Dear all, 
> I recently found the thread: Something wrong with ChartManifold when 
> chartdim=1:
> https://groups.google.com/forum/#!topic/dealii/Sfo9xKoeRpw
> a more straight answer to this thread.
>
> I took the liberty to copy David's code and modify it, so that the upper 
> face of a square follows: y(x)=0.25*sin(PI*x) + 1. 
> The function is not the same as I first proposed, but from this example it 
> is pretty clear how to proceed. 
> The result is plotted in the attached .png, along with a variation from 
> step-10 in the deal.ii tutorial.
>
> I hope this may be of any use to somebody!
>
> Juan Carlos Araújo,
> Umeå Universitet
> El martes, 9 de febrero de 2016, 17:10:51 (UTC+1), Juan Carlos Araujo 
> Cabarcas escribió:
>>
>> Dear all,
>>
>> I am working with the wave equation with variable refractive indexes 
>> within the domain.
>> I receive meshes that come from a shape optimization routine, and my task 
>> is to
>> run on arbitrarily curved elements resulting from the optimization.
>> Step-10 and others show how to use mappings and Manifolds to curve 
>> elements following basic shapes like circles, cylinders etc.
>>
>> My guess is that it should be possible in deal.II to define my own 
>> Manifolds based on how I want to bend my edges (optimization routine).
>>
>> I have prepared a very basic example on the matlab code attached that has 
>> the info to plot the attached figure. There, the edges follow the equation 
>> |x|^p + |y|^p=1, with p=6.
>>
>> It would be very illustrative to apply a mapping like what is done in 
>> step-10 on a circle, but for the presented case. I would appreciate any 
>> advise on how to achieve this, hopefully wth this example =)
>>
>> Thanks in advance!
>>
>> Juan Carlos Araújo-Cabarcas
>> Umeå Universitet.
>>
>>
>>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2001 - 2014 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Authors: Wolfgang Bangerth, Ralf Hartmann, University of Heidelberg, 2001
 * Modified from step-10 by Juan Carlos Araújo, 7/apr/2017
 */

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 

#include 
#include 
#include 

#include 

#define PI 3.14159265358979323846

namespace Step10
{
  using namespace dealii;
	
	/* ** READ ME 
		 Illustration of the use of ChartManifold, for generating a mapping 
		 of an ellipse starting from HyperBall.
		 Juan Carlos Araújo, 7/apr/2017
	*/
	
	Point<2> to_manifold (Point<2> space_point, double a, double b) {
		Point<2> R = Point<2>(space_point[0],space_point[1]), p;
		double rho, theta;
		
		theta = atan2(R[1],R[0]); // theta
if (theta < 0)
  theta += 2*PI;
		
		p = Point<2> ( a*cos(theta), b*sin(theta) );
		  
		return p;
	}
	
	// Rewritten from SphericalManifold 2D
	class EllipseManifold : public ChartManifold<2, 2, 2> {
		public:	
			EllipseManifold(double a, double b);
		
			virtual Point<2>
			pull_back(const Point<2> &space_point) const;
		
			virtual Point<2>
			push_forward(const Point<2> &chart_point) const;
		
			virtual Point<2>
			get_new_point(const Quadrature<2> &quad) const;
		
		double a, b;
		
		private:
			static Point<2> get_periodicity();	
	};	
	
	
	EllipseManifold::EllipseManifold(double a, double b):
		ChartManifold<2,2,2>(EllipseManifold::get_periodicity()), a(a),b(b)
	{}
	
	Point<2>
	EllipseMani

[deal.II] Re: Arbitrary mapping or Manifold?

[Juan Carlos Araujo Cabarcas]

Dear all, 
I recently found the thread: Something wrong with ChartManifold when 
chartdim=1:
https://groups.google.com/forum/#!topic/dealii/Sfo9xKoeRpw
a more straight answer to this thread.

I took the liberty to copy David's code and modify it, so that the upper 
face of a square follows: y(x)=0.25*sin(PI*x) + 1. 
The function is not the same as I first proposed, but from this example it 
is pretty clear how to proceed. 
The result is plotted in the attached .png, along with a variation from 
step-10 in the deal.ii tutorial.

I hope this may be of any use to somebody!

Juan Carlos Araújo,
Umeå Universitet

El martes, 9 de febrero de 2016, 17:10:51 (UTC+1), Juan Carlos Araujo 
Cabarcas escribió:
>
> Dear all,
>
> I am working with the wave equation with variable refractive indexes 
> within the domain.
> I receive meshes that come from a shape optimization routine, and my task 
> is to
> run on arbitrarily curved elements resulting from the optimization.
> Step-10 and others show how to use mappings and Manifolds to curve 
> elements following basic shapes like circles, cylinders etc.
>
> My guess is that it should be possible in deal.II to define my own 
> Manifolds based on how I want to bend my edges (optimization routine).
>
> I have prepared a very basic example on the matlab code attached that has 
> the info to plot the attached figure. There, the edges follow the equation 
> |x|^p + |y|^p=1, with p=6.
>
> It would be very illustrative to apply a mapping like what is done in 
> step-10 on a circle, but for the presented case. I would appreciate any 
> advise on how to achieve this, hopefully wth this example =)
>
> Thanks in advance!
>
> Juan Carlos Araújo-Cabarcas
> Umeå Universitet.
>
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2001 - 2014 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Authors: Wolfgang Bangerth, Ralf Hartmann, University of Heidelberg, 2001
 * Modified from step-10 by Juan Carlos Araújo, 5/apr/2017
 */

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 

#include 
#include 
#include 

#include 
#define PI 3.14159265358979323846

namespace Step10
{
  using namespace dealii;
	
	// Description of the upper face of the square
  class UpperManifold: public ChartManifold<2,2,1> {
  public:
  virtual Point<1>
  pull_back(const Point<2> &space_point) const;
  virtual Point<2>
  push_forward(const Point<1> &chart_point) const;
  };

  Point<1> UpperManifold::pull_back(const Point<2> &space_point) const {
  Point<1> p(space_point[0]);
  return p;
  }

  Point<2> UpperManifold::push_forward(const Point<1> &chart_point) const {
  double x = chart_point[0];
  // return Point<2> (x, -2*x*x + 2*x + 1);   // Parabole 
  return Point<2> ( x, 0.25*sin ( PI*x) + 1); // Trigonometric
  }

  template 
  void compute_pi_by_area ()
  {
std::cout << "Computation of Pi:" << std::endl
  << "==" << std::endl;
			
			// The case p=1, is trivial: no bending!
		  for (unsigned int degree=2; degree<11; ++degree) {
		  std::cout << "Degree = " << degree << std::endl;

		  Triangulation triangulation;

		  // new lines
	unsigned int curved_faces_label=3; // see GridGenerator::hyper_rectangle, colorize=true
	
	GridGenerator::hyper_cube (triangulation, 0.0, 1.0, true);
	
	static const UpperManifold boundary;
	triangulation.set_manifold ( curved_faces_label, boundary );
		  
		  const MappingQ mapping (degree, true);
	const QGauss quadrature(degree);
		  
	const FE_Q dummy_fe (QGaussLobatto<1>(degree + 1));
		  DoFHandler dof_handler (triangulation);

[deal.II] MappingQ, and minimal test for reproducing Benchmark

[Juan Carlos Araujo Cabarcas]

Dear all,

I'm am trying to reproduce with my implementation, the results in the 
Photonic Crystal computations performed in [1]. Here, the author uses a 
grid with an inner disk with radius R=0.475, and for FEM it is used the 
software Concepts [2] that implements curvilinear elements denoted Blending 
technique, or transfinite interpolation in quadrilaterals [3], [4,p. 144]. 
Exponential convergence was obtained for the p-FEM version.

By using the deal.II SphericalManifold, I can only get with the attached 
code up to R<0.46, before I get the exception: 

"The image of the mapping applied to cell with center [0.446317 -0.206317] 
is distorted. The cell geometry or the mapping are invalid, giving a 
non-positive volume fraction of -0.000129188 in quadrature point 54."

I am attaching a minimal test case (modified from step-10) that reproduces 
the errors:
mesh_test.cc
and the code for the same mesh used in [1], adapted to deal.II standards in:
unit_cell.h

The question is then if it is possible to get this example running for 
R=0.475, in order to reproduce the results in [1]. Is there something I am 
doing wrong that can be improved? Can this behaviour be explained by round 
off's?

Furthermore, is anyone implementing transfinite interpolation [3,4] in 
deal.II?

Thanks in advance, any help is greatly appreciated!

References:
[1] Engström, C.,Spectral approximation of quadratic operator polynomials 
arising in photonic band structure calculations, Numer. Math. (2014) 126: 
413.
[2] Concepts web page: 
http://wiki.math.ethz.ch/Concepts/WebRoot?redirectedfrom=Concepts.WebHome
[3] W. J. Gordon, C. A. Hall, Construction of curvilinear coordinate 
systems and application to mesh generation Int. J. Num. Meth. Eng., Vol. 7 
(1973), pp. 461-477
[4] Pavel Solin, Karel Segeth, Ivo Dolezel, Higher-Order Finite Element 
Methods, Studies in Advanced Mathematics, CRC Press, 2003.


-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2001 - 2014 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Authors: Wolfgang Bangerth, Ralf Hartmann, University of Heidelberg, 2001
 */

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 

#include 
#include 
#include 

// added lines
#include "unit_cell.h"
#include 

namespace Step10
{
  using namespace dealii;

  const long double pi = 3.141592653589793238462643;

  template 
  void gnuplot_output()
  {
std::cout << "Output of grids into gnuplot files:" << std::endl
  << "===" << std::endl;

Triangulation triangulation;
GridGenerator::hyper_ball (triangulation);
static const SphericalManifold boundary;
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold (0, boundary);

for (unsigned int refinement=0; refinement<2;
 ++refinement, triangulation.refine_global(1))
  {
std::cout << "Refinement level: " << refinement << std::endl;

std::string filename_base = "ball";
filename_base += '0'+refinement;

for (unsigned int degree=1; degree<4; ++degree)
  {
std::cout << "Degree = " << degree << std::endl;

const MappingQ mapping (degree);

GridOut grid_out;
GridOutFlags::Gnuplot gnuplot_flags(false, 30);
grid_out.set_flags(gnuplot_flags);

std::string filename = filename_base+"_mapping_q";
filename += ('0'+degree);
filename += ".dat";
std::ofstream gnuplot_file (filename.c_str());

grid_out.write_gnuplot (triangulation, gnuplot_file, &mapping);
  }
std::cout << std::endl;
  }
  }

  template 
  void compute_pi_by_area ()
  {
std::cout << "Computation of Pi by the area:" << std::endl
  << "==" << std::endl;

// const QGauss quadrature(4);
			
		for (double R=0.45; R < 0.5; R += 0.005)
		  for (unsigned int degree=1; degree

[deal.II] Re: periodic boundary conditions

[Juan Carlos Araujo Cabarcas]

Hej, I found the issue ... I forgot the line

 constraints.distribute ( sol );

after the solution routine ... so it was just mainly a visualization issue. 
Thanks for the quick response!

El jueves, 22 de septiembre de 2016, 19:26:28 (UTC+2), Daniel Arndt 
escribió:
>
> Juan,
>
> [...]
>> 1) My first question is about how the "direction" works, is it 1 
>> correspond to "x" and 2 corresponds to "y"?
>>
> No, direction 0/1/2 means that the support_points of the dofs to be 
> identified just differ in the x-/y-/z-component.
> What you are doing should be perfectly fine.
>  
>
>> 2) When I do this I obtain fairly reasonable solutions, however I am 
>> getting a weird Dirichlet-type of boundary conditions on boundaries 1 and 3 
>> (see the attached figures), that distort the fields.
>> I seem unable to find what is wrong, so I wonder if anyone has 
>> experienced this kind of behavior.
>>
> No, I haven't observed this problem. Can you check in your 
> ConstraintMatrix that you don't use any other constraints? For doing so use 
> ConstraintMatrix.print(std::cout)
> and check if there is a DoF that is not constrained to another one.
>  
> Best,
> Daniel
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Fast FE evaluation over several points and several functions

[Juan Carlos Araujo Cabarcas]

Thanks for the fast reply, this is exactly what I needed!

El miércoles, 22 de junio de 2016, 19:38:35 (UTC+2), bangerth escribió:
>
>
> > I am solving an eigenvalue problem similar than step-36. After solving 
> for the 
> > eigenpairs, I evaluate the eigenfunctions in the standard way: 
> > 
> >  VectorTools::point_value ( mapping, dof_handler, efun[m], 
> > q_points[j], Uq ); 
> > 
> > where: efun[m] is the m-th eigenfunction from step-36, 
> > q_points[j] are selected quadrature points, 
> > Uq is where I store the FE evaluation. 
> > 
> > As expected, this operation is awfully slow! it takes seconds for a 
> single 
> > point evaluation with a decent discretization and having m 
> eigenfunctions 
> > makes it worse! 
>
> Yes. But there is a much more efficient way to do this: 
>FEValues fe_values (...); 
>std::vector sol_at_q_points (...); 
>for (cell=...) 
>  { 
> fe_values.reinit (cell); 
> fe_values.get_function_values (efun[m], sol_at_q_points); 
>
> This gives you the values of efun[m] at all quadrature points on the 
> current 
> cell at once, and this approach is efficient and independent of the 
> overall 
> size of the problem. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth   email:bang...@math.tamu.edu 
>  
>  www: http://www.math.tamu.edu/~bangerth/ 
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Fast FE evaluation over several points and several functions

[Juan Carlos Araujo Cabarcas]

Dear all, 

I am solving an eigenvalue problem similar than step-36. After solving for 
the eigenpairs, I evaluate the eigenfunctions in the standard way:

VectorTools::point_value ( mapping, dof_handler, efun[m], 
q_points[j], Uq );

where: efun[m] is the m-th eigenfunction from step-36, 
q_points[j] are selected quadrature points, 
Uq is where I store the FE evaluation.

As expected, this operation is awfully slow! it takes seconds for a single 
point evaluation with a decent discretization and having m eigenfunctions 
makes it worse!

1) I wonder if there is a way to evaluate the whole vector q_points with a 
single (and clever) call instead of looping on j and calling 
VectorTools::point_value(... , q_points[j], ...).

I remember when coding basic FEM in matlab, I had loops over several points 
in order to reduce the overhead of evaluating the FE

for c in cells
for x in q_points
if x is in cell
evaluate FE: evaluate FE:  loop over shape functions with 
support in c, 
U(x)=sum ...
end
end
end

2) My issue today goes beyond evaluating in several points ... I also 
require to evaluate several FE vectors (eigenfunctions). Intuitively, one 
would evaluate the FE like:

for c in cells
for x in q_points
if x is in cell
for m ... in eigenfunctions
evaluate FE:  loop over shape functions with support in c, 
Um(x)=sum ...
end
end
end
end

Any ideas how to achieve this?
Thanks in advance.

Juan Carlos Araújo
Umeå Universitet


-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: precompiling independent classes as a library and link it to the main project

[Juan Carlos Araujo Cabarcas]

Thanks for the quick response ... I will give it a try and come back to you.

El martes, 31 de mayo de 2016, 5:38:37 (UTC+2), David Wells escribió:
>
> Hi Juan,
>
> I use this approach (shared object library for pieces, executable for the 
> solver) approach frequently and it is a bit tricky to set up.
>
> What I usually do is have three different CMakeLists.txt files: one to 
> control everything, one for the library, and one for the solver. I 
> implemented an example of this in the cdr code gallery entry:
>
> https://github.com/dealii/code-gallery/tree/master/cdr
>
> Let me know if this helps.
>
> Thanks,
> David Wells
>
> On Monday, May 30, 2016 at 10:35:14 AM UTC-4, Juan Carlos Araujo Cabarcas 
> wrote:
>>
>> Dear all, 
>>
>> I have several independent classes that use deal.II objects, and I would 
>> like to pre-compile them as shared libraries to be linked to the main 
>> executable.
>> The idea is to be able to do something similar than what is described 
>> here:
>> http://www.cprogramming.com/tutorial/shared-libraries-linux-gcc.html
>> but as I am not really used to the cmake environment I have not been able 
>> to reproduce a similar example with deal.ii.
>>
>> As a test example I would like to remove the class BoundaryValues from 
>> step-4:
>> https://dealii.org/8.4.0/doxygen/deal.II/step_4.html
>> and save it as a file bval.cc pre-compile it to generate a .so file to be 
>> linked in the main executable.
>>
>> Any ideas how to achieve this?
>>
>> Thanks in advance,
>>
>> Juan Carlos Araújo
>> PhD student at Umeå Universitet
>>
>>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] precompiling independent classes as a library and link it to the main project

[Juan Carlos Araujo Cabarcas]

Dear all, 

I have several independent classes that use deal.II objects, and I would 
like to pre-compile them as shared libraries to be linked to the main 
executable.
The idea is to be able to do something similar than what is described here:
http://www.cprogramming.com/tutorial/shared-libraries-linux-gcc.html
but as I am not really used to the cmake environment I have not been able 
to reproduce a similar example with deal.ii.

As a test example I would like to remove the class BoundaryValues from 
step-4:
https://dealii.org/8.4.0/doxygen/deal.II/step_4.html
and save it as a file bval.cc pre-compile it to generate a .so file to be 
linked in the main executable.

Any ideas how to achieve this?

Thanks in advance,

Juan Carlos Araújo
PhD student at Umeå Universitet

-- 
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.
For more options, visit https://groups.google.com/d/optout.