[deal.II] FE_Q(p) with p>1: Global coordinates of midnodes

[Simon]

Dear All,

I am implementing a local postprocessing SPR-Approach for a mechanical 
problem, in which I sloppy speaking need the global coordinates of all 
nodes.

I use Elements of type FE_Q(p). 
For p=1 there is no problem: I can simply loop over all cells and vertices 
and use cell->vertex(...) in order to get the global coordinates of 
vertex(...). 
Being p=2 for instance, this approach does not work anymore. In this case 
is there an easy way to figure out the coordinates of midnodes as well?

Thanks for help!

Best
Simon

-- 
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/3d87b2c4-61a4-4125-982d-f7341f44492an%40googlegroups.com.

Re: [deal.II] FE_Q(p) with p>1: Global coordinates of midnodes

[Simon]

Hi Paras,

everything worked fine, thanks a lot!

Best 
Simon

Paras Kumar schrieb am Montag, 19. April 2021 um 12:28:37 UTC+2:

> Hi Simon,
>
> This could be helpful: 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
>
> Here's a code snippet for Option-1, which I usually use.
>
> const auto &   feSystem   = dofHandler.get_fe();
> dealii::MappingQ cellMapping(feSystem.degree);
> const auto &   referenceCellSupportPoints =  
> feSystem.get_unit_support_points();
> std::vector> 
> realCellSupportPoints(referenceCellSupportPoints.size());
> unsigned int supportPointCtr = 0;
> for (const auto &supportPoint : referenceCellSupportPoints)
>   {
> realCellSupportPoints[supportPointCtr] =
>   cellMapping.transform_unit_to_real_cell(cellIter, supportPoint);
> ++supportPointCtr;
>   }
>
>
> Best regards,
> Paras
>
> On Mon, Apr 19, 2021 at 12:14 PM Simon  wrote:
>
>> Dear All,
>>
>> I am implementing a local postprocessing SPR-Approach for a mechanical 
>> problem, in which I sloppy speaking need the global coordinates of all 
>> nodes.
>>
>> I use Elements of type FE_Q(p). 
>> For p=1 there is no problem: I can simply loop over all cells and 
>> vertices and use cell->vertex(...) in order to get the global coordinates 
>> of vertex(...). 
>> Being p=2 for instance, this approach does not work anymore. In this case 
>> is there an easy way to figure out the coordinates of midnodes as well?
>>
>> Thanks for help!
>>
>> Best
>> Simon
>>
>> -- 
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see 
>> https://groups.google.com/d/forum/dealii?hl=en
>> --- 
>> You received this message because you are subscribed to the Google Groups 
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an 
>> email to dealii+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/3d87b2c4-61a4-4125-982d-f7341f44492an%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/3d87b2c4-61a4-4125-982d-f7341f44492an%40googlegroups.com?utm_medium=email&utm_source=footer>
>> .
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/84f26166-9200-4d0e-aea5-580c15933f3dn%40googlegroups.com.

[deal.II] problem with CylindricalManifold<3> after removing cells from a triangulation

[Simon]

Dear all,
 
I created a mesh using GridGenerator::hyper_cube_with_cylindrical_hole and 
removed after that half of the total number of the cells in order to make 
use of symmetry. 
Then I set a SphericalManifold<2> (2D) respectively CylindricalManifold<3> 
(3D) in order to get a true circle (2D) respectively a true cylinder (3D) 
after global and/or adaptive refining. 

Without removing some cells, i.e. just using the 
GridGenerator::hyper_cube_with_cylindrical_hole without using symmetry, 
both 2D and 3D variants work fine. 
With removing some cells, i.e. using symmetry, only in 2D I get a true 
circle as result. In 3D however, the SphericalManifold<3> does not work 
anymore. (I can adaptively refine the inner hole, but the result is not a 
true cylinder anymore which was the case without removing the cells).

In the following is a snippet of my code. It´s hard for me to find out what 
the problem actually is, since in 2D the manifold description worked after 
removing some cells. 

Thanks for helping!

Best
Simon




const double outer_radius = 1.0;
const double inner_radius = 0.5;
const Point center;

GridGenerator::hyper_cube_with_cylindrical_hole(triangulation_tmp,
inner_radius,
 outer_radius,
0.5,
1,
 false 
/*boundary_id_inner_hole is set to 1*/);

std::set::active_cell_iterator> cells_to_remove;
typename Triangulation::active_cell_iterator cell 
=triangulation_tmp.begin_active(),
endc = triangulation_tmp.end();

for(; cell!=endc; cell++)
 {
  if(cell->center()[1]<0)
   {
cells_to_remove.insert(cell);
}
  }
GridGenerator::create_triangulation_with_removed_cells(triangulation_tmp, 
cells_to_remove, triangulation);

std::unique_ptr> ptr_manifold=nullptr;

if(dim==2)
{
ptr_manifold = std::make_unique>(center);
}
else if(dim==3)
{
ptr_manifold = std::make_unique>(dim-1);
}
else
{
throw std::runtime_error("only allowed for dim == 2 or dim == 3");
}
types::boundary_id boundary_id_inner_hole=1;
types::manifold_id manifold_id_inner_hole=1;
   
triangulation.set_all_manifold_ids_on_boundary(boundary_id_inner_hole,manifold_id_inner_hole);

triangulation.set_manifold (manifold_id_inner_hole, *ptr_manifold);

..

-- 
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/ab270946-bbe3-4952-a862-b7fd2e6ccd3bn%40googlegroups.com.

[deal.II] Re: problem with CylindricalManifold<3> after removing cells from a triangulation

[Simon]

By the way, I forgot to precise my question in more detail:

In the code snippet above I set the corresponding Manifold Descriptions for 
dim=2 (spherical) and dim=3 (cylindrical) on the inner boundary created by  
*GridGenerator::hyper_cube_with_cylindrical_hole* and removed after that 
some cells by calling 
*GridGenerator::create_triangulation_with_removed_cells* . 
However, for the 3D-Case, after adaptively refining the inner hole there is 
still a FlatManifold attached to it and not a CylindricalManifold. For 
dim=2, I get the desired SphericalManifold and everything is fine. 
I am pretty sure that I set the manifold_id on the inner boundary 
correctly, because after visualizing the mesh the cells are actually 
adaptively refined around the inner boundary, but as I said a FlatManifold 
is used in the 3D-Case. 

-So does *GridGenerator::hyper_cube_with_cylindrical_hole* attach some 
important "information" in terms of id_s,... to my triangulation which my 
new triangulation created by  
*GridGenerator::create_triangulation_with_removed_cell* does not have?
-Or does the CylindricalManifold only work for closed boundaries? Because 
without calling *GridGenerator::create_triangulation_with_removed_cells *also 
in the 3D-Case I get the desired CylindricalManifold.

Best
Simon
Simon schrieb am Freitag, 23. April 2021 um 22:02:41 UTC+2:

> Dear all,
>  
> I created a mesh using GridGenerator::hyper_cube_with_cylindrical_hole and 
> removed after that half of the total number of the cells in order to make 
> use of symmetry. 
> Then I set a SphericalManifold<2> (2D) respectively CylindricalManifold<3> 
> (3D) in order to get a true circle (2D) respectively a true cylinder (3D) 
> after global and/or adaptive refining. 
>
> Without removing some cells, i.e. just using the 
> GridGenerator::hyper_cube_with_cylindrical_hole without using symmetry, 
> both 2D and 3D variants work fine. 
> With removing some cells, i.e. using symmetry, only in 2D I get a true 
> circle as result. In 3D however, the SphericalManifold<3> does not work 
> anymore. (I can adaptively refine the inner hole, but the result is not a 
> true cylinder anymore which was the case without removing the cells).
>
> In the following is a snippet of my code. It´s hard for me to find out 
> what the problem actually is, since in 2D the manifold description worked 
> after removing some cells. 
>
> Thanks for helping!
>
> Best
> Simon
>
>
>
>
> const double outer_radius = 1.0;
> const double inner_radius = 0.5;
> const Point center;
>
> GridGenerator::hyper_cube_with_cylindrical_hole(triangulation_tmp,
> inner_radius,
>  outer_radius,
> 0.5,
> 1,
>  false 
> /*boundary_id_inner_hole is set to 1*/);
>
> std::set::active_cell_iterator> 
> cells_to_remove;
> typename Triangulation::active_cell_iterator cell 
> =triangulation_tmp.begin_active(),
> endc = triangulation_tmp.end();
>
> for(; cell!=endc; cell++)
>  {
>   if(cell->center()[1]<0)
>{
> cells_to_remove.insert(cell);
> }
>   }
> GridGenerator::create_triangulation_with_removed_cells(triangulation_tmp, 
> cells_to_remove, triangulation);
>
> std::unique_ptr> ptr_manifold=nullptr;
> 
> if(dim==2)
> {
> ptr_manifold = std::make_unique>(center);
> }
> else if(dim==3)
> {
> ptr_manifold = std::make_unique>(dim-1);
> }
> else
> {
> throw std::runtime_error("only allowed for dim == 2 or dim == 3");
> }
> types::boundary_id boundary_id_inner_hole=1;
> types::manifold_id manifold_id_inner_hole=1;
>
> triangulation.set_all_manifold_ids_on_boundary(boundary_id_inner_hole,manifold_id_inner_hole);
>
> triangulation.set_manifold (manifold_id_inner_hole, *ptr_manifold);
> 
> ..
>
>

-- 
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/d326f1ec-97db-419c-a055-f727ebe55454n%40googlegroups.com.

Re: [deal.II] problem with CylindricalManifold<3> after removing cells from a triangulation

[Simon]

Hello Luca,

your suggestion fixed my problem, so thank you very much!

Just to make sure that I understand it correctly: 
In my code snippet from above I call the member function 
"triangulation.set_all_manifold_ids_on_boundary(1,1). 
I actually thought that this call does the same thing as you suggested, 
i.e. cell->face(f)->set_all_manifold_ids(1).
But I guess the problem is that in my version I only set the manifold ids 
of the face and not its four edges as well, is it?

Best
Simon



luca@gmail.com schrieb am Montag, 26. April 2021 um 09:56:17 UTC+2:

> There may be a problem with the way boundary ids are set.
>
> Can you try the following?
>
> after creating the grid with removed cells, loop over all cells and all 
> faces, and if at boundary with boundary id == 1, then call
>
> cell->face(f)->set_all_manifold_ids(1);
>
> notice the “_all_”, i.e., it is not “set_manifold_id(1)”. This ensures you 
> set manifold ids also on edges.
>
> It may be that boundary ids are not set correctly on edges when creating 
> the triangulation with removed cells. In this case this is a bug (simple to 
> fix, btw), that you only see in 3d, because refinement starts from edges.
>
> L.
>
>
>
> > On 26 Apr 2021, at 9:36, Simon  wrote:
> > 
> > By the way, I forgot to precise my question in more detail:
> > 
> > In the code snippet above I set the corresponding Manifold Descriptions 
> for dim=2 (spherical) and dim=3 (cylindrical) on the inner boundary created 
> by GridGenerator::hyper_cube_with_cylindrical_hole and removed after that 
> some cells by calling 
> GridGenerator::create_triangulation_with_removed_cells . 
> > However, for the 3D-Case, after adaptively refining the inner hole there 
> is still a FlatManifold attached to it and not a CylindricalManifold. For 
> dim=2, I get the desired SphericalManifold and everything is fine. 
> > I am pretty sure that I set the manifold_id on the inner boundary 
> correctly, because after visualizing the mesh the cells are actually 
> adaptively refined around the inner boundary, but as I said a FlatManifold 
> is used in the 3D-Case. 
> > 
> > -So does GridGenerator::hyper_cube_with_cylindrical_hole attach some 
> important "information" in terms of id_s,... to my triangulation which my 
> new triangulation created by 
> GridGenerator::create_triangulation_with_removed_cell does not have?
> > -Or does the CylindricalManifold only work for closed boundaries? 
> Because without calling 
> GridGenerator::create_triangulation_with_removed_cells also in the 3D-Case 
> I get the desired CylindricalManifold.
> > 
> > Best
> > Simon
> > Simon schrieb am Freitag, 23. April 2021 um 22:02:41 UTC+2:
> > Dear all,
> > 
> > I created a mesh using GridGenerator::hyper_cube_with_cylindrical_hole 
> and removed after that half of the total number of the cells in order to 
> make use of symmetry. 
> > Then I set a SphericalManifold<2> (2D) respectively 
> CylindricalManifold<3> (3D) in order to get a true circle (2D) respectively 
> a true cylinder (3D) after global and/or adaptive refining. 
> > 
> > Without removing some cells, i.e. just using the 
> GridGenerator::hyper_cube_with_cylindrical_hole without using symmetry, 
> both 2D and 3D variants work fine. 
> > With removing some cells, i.e. using symmetry, only in 2D I get a true 
> circle as result. In 3D however, the SphericalManifold<3> does not work 
> anymore. (I can adaptively refine the inner hole, but the result is not a 
> true cylinder anymore which was the case without removing the cells).
> > 
> > In the following is a snippet of my code. It´s hard for me to find out 
> what the problem actually is, since in 2D the manifold description worked 
> after removing some cells. 
> > 
> > Thanks for helping!
> > 
> > Best
> > Simon
> > 
> > 
> > 
> > 
> > const double outer_radius = 1.0;
> > const double inner_radius = 0.5;
> > const Point center;
> > 
> > GridGenerator::hyper_cube_with_cylindrical_hole(triangulation_tmp,
> > inner_radius,
> > outer_radius,
> > 0.5,
> > 1,
> > false /*boundary_id_inner_hole is set to 1*/);
> > 
> > std::set::active_cell_iterator> 
> cells_to_remove;
> > typename Triangulation::active_cell_iterator cell 
> =triangulation_tmp.begin_active(),
> > endc = triangulation_tmp.end();
> > 
> > for(; cell!=endc; cell++)
> > {
> > if(cell->center()[1]<0)
> > {
> > cells_to_remove.insert(cell);
> > }
> > }
> > 
> GridGenerator::create_triangulation_with_removed_cells(triangulation_tm

[deal.II] GridTools::find_cells_adjacent_to_vertex: how to avoid multiple calling ?

[Simon]

Dear all,

the mentioned function is exactly doing the right thing for me, i.e. 
returning an vector of iterators to all cells, which share this vertex. 
After calling this function I will loop over these cells and compute 
certain things ...

My problem is that I have to work with two different dof_handlers and also 
need some quantities from FEValues of both DoF_Handlers on my neighbor 
cells. So far I handled this like this:

auto cells_ref = GridTools::find_cells_adjacent_to_vertex(dof_handler_ref, 
/*vertex*/);
auto cells_tmp = GridTools::find_cells_adjacent_to_vertex(dof_handler_tmp, 
/*vertex*/); 

I am not sure how "expensive" this function actually is. I have to call 
this function for all vertices of my triangulation, so currently I call it 
"2*n_vertices()" times.

I´d like to call this function only once for a given vertex, i.e. 
"n_vertices()"times.
Is there a way to realize this in combination with two different 
dof_handlers?

Thanks for helping!

Best
Simon


-- 
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/4de4a651-1092-47cc-9a72-b077a9e7169an%40googlegroups.com.

[deal.II] calculation of the L2-norm using "integrate_difference" - mismatch with own test case

[Simon]

Dear all,

I want to compute the L2-norm for a large number of scalar 
history-variables. The solution vectors for all of them are stored in the 
following container: 
"std::vector>(num_of_history_variables)".
So far I called integrate_difference " num_of_history_variables" times, 
which is quite time consuming. Thus, I´d like to compute the L2-norm by 
myself to get a second time measurement and may save in the first call some 
time reducing information which integrate_differnce doesn´t provide for 
subsequent calls.

In order to do so, I compared for a arbitrary history variable the vector 
of cell errors, calculated from integrate_differnce, with those cell errors 
of my own implementation. Unfortunately I don´t come up with the same cell 
error(s). 

In the following is a snippet of my code, some remarks for a better 
understanding: 
(I access my history variables via a CellDataStorage object. The Variable 
"comp_0" is the value at the QP of my current mesh, the variable 
"comp_0_ref" is the corresponding value of the reference solution on the 
same QP).

Functions::FEFieldFunction fe_field(dof_handler_scalar_ar, 
reference_solution_n[0]);

cell=dof_handler_scalar.begin_active();
for(; cell!=endc; cell++)
{
fe_values.reinit(cell);

std::vector>> lqph =  
quadrature_point_history.get_data(cell);

double error_cell=0;

 for(unsigned int q=0; q projection_variables = 
lqph[q]->create_vector_for_projection();
 double comp_0 = projection_variables(0);

 error_cell += ((comp_0 - comp_0_ref)*(comp_0 - comp_0_ref)*JxW);
 }
std::cout<<"Local cell error: "

[deal.II] Re: calculation of the L2-norm using "integrate_difference" - mismatch with own test case

[Simon]

By the way, I forgot to mention how I acutally compute the variable 
"comp_0_ref", which is potentially the place where the difference may come 
from.
I do this by simply asking the initialized FE Field function: 
double comp_0_ref = fe_field.value(fe_values.quadrature_point(q));

For my understanding of the L2-Norm this is acutally the value that I need, 
i.e. I need the value of the reference solution at the qps of my current 
solution, do I?

-Simon

Simon schrieb am Mittwoch, 12. Mai 2021 um 16:47:42 UTC+2:

>
> Dear all,
>
> I want to compute the L2-norm for a large number of scalar 
> history-variables. The solution vectors for all of them are stored in the 
> following container: 
> "std::vector>(num_of_history_variables)".
> So far I called integrate_difference " num_of_history_variables" times, 
> which is quite time consuming. Thus, I´d like to compute the L2-norm by 
> myself to get a second time measurement and may save in the first call some 
> time reducing information which integrate_differnce doesn´t provide for 
> subsequent calls.
>
> In order to do so, I compared for a arbitrary history variable the vector 
> of cell errors, calculated from integrate_differnce, with those cell errors 
> of my own implementation. Unfortunately I don´t come up with the same cell 
> error(s). 
>
> In the following is a snippet of my code, some remarks for a better 
> understanding: 
> (I access my history variables via a CellDataStorage object. The Variable 
> "comp_0" is the value at the QP of my current mesh, the variable 
> "comp_0_ref" is the corresponding value of the reference solution on the 
> same QP).
>
> 
> Functions::FEFieldFunction fe_field(dof_handler_scalar_ar, 
> reference_solution_n[0]);
>
> cell=dof_handler_scalar.begin_active();
> for(; cell!=endc; cell++)
> {
> fe_values.reinit(cell);
>
> std::vector>> lqph =  
> quadrature_point_history.get_data(cell);
>
> double error_cell=0;
> 
>  for(unsigned int q=0; q  {
>  const double JxW = fe_values.JxW(q);
>
>  Vector projection_variables = 
> lqph[q]->create_vector_for_projection();
>  double comp_0 = projection_variables(0);
>
>  error_cell += ((comp_0 - comp_0_ref)*(comp_0 - comp_0_ref)*JxW);
>  }
> std::cout<<"Local cell error: "< }
>
> --
>
> I had a look in the implementation of the integrate_difference function 
> with case == L2 as well.
> However I could not identify any differences to my own implementation.
>
> My question is if there is a mistake in my code snippet, i.e in the way 
> how I calculate the integrals for the cell errors? 
> On each cell, I basically multiply the squared abs-value of the difference 
> with the JxW-Value, sum those values up and finally take the square root.
>
> Thanks in advance for helping!
>
> Best
> Simon
>

-- 
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/946ac5ef-8ce9-407e-97f0-f39297b8481en%40googlegroups.com.

[deal.II] projectin of gradients of the discrete solution to nodes - two different discontinuous approaches

[Simon]

Dear all,

I am projecting gradients of my displacement solution, for instance 
stresses, which are only available to me at the quadrature points, to the 
nodes in order to get fields. (However I project many more variables, so I 
transformed my problem in projecting many scalar variables, but anyway this 
is not important regarding my question).

-I am doing the projection once by solving a global minimization problem, 
i.e. minimizing the integral over the whole domain of the squared 
difference diff = "values at qps - values at nodes". This leads to the 
typical mass matrix (phi_i, phi_j) which is also provided in the 
VectorTools namespace.
I solve this min.-problem once with a continuous FE_Q und once with a 
discontinuous FE_DGQ, i.e. I call the very same function body, I just give 
different DoFHandlers as input.

-A second approach for the same discontinuous FE_DGQ is intruduced in 
step-18 (extensions) by using 
compute_projection_from_quadrature_points_matrix(), i.e. a local method. 
(My poly_degree is constructed in a way that I have always as many 
dofs_per_cell as qps.)

My question is if the output, i.e. the nodal values, of both discontinuous 
approaches 
(1. given a FE_DGQ to the global min.-problem vs. 2. the local 
compute_projection_from_quadrature_points_matrix() approach)
should be the same, except of numerical errors?

I am aware that I solve a linear system in the first case whereas in the 
second case not, but my idea was the following: The first approach 
minimizes the squared difference, but if I have as many dofs_per_cell as 
qps then of course the squared difference can be minimized to zero for each 
qp, i.e. the qp values can directly be transferred to the nodes. And if my 
understanding is correct, this is exactly what the second approach does. So 
in the first approach I do not do it "directly" but due to its definition 
and the number of dofs it should do the same as second.
Of course I compared this approaches in my program. However depending on 
the mesh size there are deviations up to percent, with finer meshes this 
difference reduces. 
Can this deviation be argumented away with the standard argument "this is 
the numerics..." or is there is a mathematical difference and both 
approaches do something different?

Thanks for the input!

Best
Simon

-- 
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/f04edf1d-cc06-4732-bb3d-d3e74162ed98n%40googlegroups.com.

[deal.II] SparseDirectUMFPACK - member function "initialize" takes very long

[Simon]

Dear all,

I use the SparseDirectUMFPACK solver for doing some post-processing steps.
I also measure some cpu_times in my program in order to compare different 
methods regarding efficiency. For my question only the following three 
lines of code from my program are relevant:
---
timer.enter_section("Initialize Direct Solver");
sparse_direct_solver.initialize(global_matrix);
timer.leave_subsection("Initialize Direct Solver");
---
So I only measure the time which is neeed for executing the member function 
"initialize":
->Solving a scalar problem with *dim=2* and *263425 DoFs* the above call 
takes about 1.3 seconds, whereas solving the linear system takes only 0.8 
seconds.
->Solving the same scalar problem with *dim=3* and only *72369 DoFs* the 
above call takes about 4.4 seconds and solving the linear system about 3.2 
seconds.

-My first question is what this function actually does since it takes 
longer than assembling and solving the system together? 
The documentation says: "This function does nothing. It is only here to 
provide a interface consistent with other sparse direct solvers."
By having a look in the function body I´ve seen that this function acutally 
does nothing, i.e. the function body is empty.

-My second question is why for my dim=3 case the function call takes much 
longer as in the dim=2 case? In the latter I have much more DoFs.
My guess was that the number of DoFs is the essential quantity but this is 
obviously not the case.

Any input would be greatly appreciated!

Best
Simon

-- 
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/c23feb6d-faf9-4b52-a25d-f639d309ba13n%40googlegroups.com.

Re: [deal.II] Re:

[Simon]

Yes, actually both schemes converge to the same solution. 
Maybe I should do some postprocessing with a visualiuation in order to 
figure out what the lumped integration does with the values at the qps.

I´ve got one more basic question ragarding the determination of the 
convergence rate as it is done in many tutorials (step 7, step 51,...). For 
that purpose I watched your new videos (3.9, 3.91, 3.93, 3.95) where you 
derived the error estimation for the laplace equation in certain norms. You 
mentioned in there that for many PDEs one comes up with the same results, 
i.e. that the finite element error is at least as good as the interpolation 
error *times some equation-dependent constant*. 
I want to determine the convergence rate for the elasticity equations with 
the PDE: div(stress-tensor) = 0. 
-Since I determine the convergence rate, the before mentioned *constant* in 
the inequality cancels out, i.e. I should get the same convergence rate as 
for the laplace equation. Is that right?

-I actually solve the nonlinear elasticity equations (hyperelasticity); The 
finite element spaces, test functions,... are the same as for linear 
elasticity, but the stresses depend now nonlinear on the gradient of the 
displacements. My question is if the nonlinearity changes the error 
estimates?
Let´s assume I determine the convergence rate in the L2-norm for linear 
elasticity, which is in the best case 2 for for p=1. Should I also get 2 
for the nonlinear elasticity equations?

I know that these questions are not directly related to dealii, but I guess 
you have a deep knowledge in that area. I would appreciate if you helped my 
answering that.

Best
Simon

Wolfgang Bangerth schrieb am Sonntag, 6. Juni 2021 um 01:58:39 UTC+2:

> On 6/4/21 9:14 AM, Simon Wiesheier wrote:
> > 
> > So intuitively speaking, what is the effect of the lumped integration? 
> From a 
> > mathematical viewpoint there are no more couplings between the DoFs. But 
> does 
> > this lead to lower nodal values, higher nodal values or is this context 
> dependent?
>
> I have to admit that I have no intuition.
>
> Do both schemes converge? To the same solution?
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/30a00e59-625d-403e-812f-cc051802deeen%40googlegroups.com.

Re: [deal.II] Re: SparseDirectUMFPACK - member function "initialize" takes very long

[Simon]

 One more question came up when I watched your video "What solver to use". 
In there you mentioned that direct solvers in 2d have a complexity of 
O(N^2), where N is the number of unknowns. 
There is an approximation for N, i.e. 
N \approx p^d *|Omega| / h^d.
So from that can one say that solving a linear system with UMFPACK has a 
complexity of O(1/h^{4}) for a fixed degree p and d=2?
I ask that because I measured the times for solving two large linear 
systems in 2d, one with 1 Million DoFs and the other one with 2 Million 
DoFs (reducing h to h/2). The former took 11 seconds, the latter 110 
seconds, i.e. they are related by a factor of 10 and not 16.

Backcalculating this, the factor of 10 would follow from a complexity of 
O(N^5/3}). 
So is this is a plausible rate or is my linear system still not big enough 
to see the rate of 16? Is the renumbering scheme, which does UMFPACK 
internally, may the source why I see a lower rate?

Best
Simon
Wolfgang Bangerth schrieb am Mittwoch, 2. Juni 2021 um 18:21:56 UTC+2:

> On 6/2/21 10:16 AM, Simon Wiesheier wrote:
> > One more question to the renumbering schemes: I called 
> > DoFRenumbering::Cuthill_McKee(dof_handler) before calling the initialize 
> > function of the direct solver.
> > 
> > So this isn´t really necessary when using UMFPACK because internally a 
> > different renumbering scheme is called anyway?
>
> Correct.
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/168e3b0a-5860-4c0a-9d70-4639d13e63a1n%40googlegroups.com.

[deal.II] writing values in the solution vector for constrained dofs

[Simon]

Dear all,

I am dealing with adaptivity, so hanging nodes are present in my mesh which 
are handled by a constraints object like in the following:
AffineConstraints...constraints;
make_hanging_node_constraints...();
constraints.close();

My solution vector is filled via a local approach. Basically I loop over 
all active vertices, use cell->vertex_dof_index() to access the dofs living 
on them and write appropriate values in the global solution vector. So I 
don't use the before mentioned constraints object to make a local-global 
distribution or something similar as it is done in most tuturial programs. 
After my loop is finished I call constraints.distribute().

So by looping over all active vertices I also visit the constrained hanging 
nodes and write some values for these dofs in the solution vector.

My question is if I can just skip the hanging nodes in my loop because the 
call "constraints.distribute()" will overwrite these values anyway?
This is at least my understanding of what constraints.distribute() does, 
i.e. computing constraind dof values from unconstrained ones. But I am not 
quite sure if problems can come up by modifyiing the constrained dofs 
beforehand. 

Best
Simon

-- 
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/d299ea70-e704-4edc-ada1-b1ebc554828bn%40googlegroups.com.

[deal.II] std::set> : missing operator(s)

[Simon]

Dear all,

I am trying to create a set of Points, i.e std::set>. 
I get a quite extensive error message when trying to insert a point in the 
set but for my understanding of the message, the only problem is the 
missing operator< for two points. 
(Well I could make use of a std::vector instead but a std::set fits much 
better for my issue since (i) I need to insert points at arbitrary 
positions quite often and second (ii) wanna make sure that a given Point is 
only once in the container.)

But irrespective of what container I use, I would like to supplement the 
Point Class by adding the 'operator<' like in the following:

template
bool Point::operator< (const Point& p)
{
//...some Asserts
return this.norm()https://www.dealii.org/current/doxygen/deal.II/base_2point_8h_source.html>> 
by simply adding the missing operator< declaration and the corresponding 
definition. So I did this but without effect, probably because all the 
files of the dealii-library are already compiled and executing my program 
via 'make run' did not see my changes in the Point Class at all. (I also 
added a dummy command, std::cout<<"...Test..."

[deal.II] solution of linear system using LAPACKFullMatrix - problems with singular matrix

[Simon]

Dear all,

I have to deal with a large number of very small linear systems A*a=b; I 
need to solve this system for each vertex of my triangulation. 
The matrix A is purely determinand by dyadic products of certain quadrature 
points.
Example: Point qp =[x,y]; Vector qp_ = [1, x, y, xy, x^2, 
y^2]^T
-> A = A + qp_*qp_^T
The number of qps is restricted to the qps of those cells which are in the 
patch. For instance there are 16 qps when using Q2 elements in 2d with 
reduced integration (QGauss(2)).

First I decided to actually invert the matrix A, 
i.e. A.invert(A), 
but the computation of the inverse introduced too much numerical errors. 
Hence I decided to switch to LAPACKFullMatrix, computed LU factorization 
and finally solved the system via the .solve() member functions. 

Switching to LAPACK made my results much better but I still reach a point 
(fine meshes) where the LU decomposition fails. For nearly each vertex the 
reciprocal condition number reaches the machine accuracy.
I guess the problem is that fine meshes lead to large differences in the x- 
and y-coordinates of  the qps (Point dummy = [x,y] = [10.0 ,0.01]) so 
that some entries in the A matrix are close to zero and others in turn 
become very large. IMHO the LU decomposition becomes more demanding if the 
entries in A differ by several orders. Here is a typical example for A 
(lowest order e-01 , highest order e+02):

1.600e+01  6.099e+01  3.798e+00  2.327e+02  1.203e+00  1.449e+01 
 6.099e+01  2.327e+02  1.449e+01  8.884e+02  4.594e+00  5.532e+01 
 3.798e+00  1.449e+01  1.203e+00  5.532e+01  4.292e-01  4.594e+00 
 2.327e+02  8.884e+02  5.532e+01  3.395e+03  1.755e+01  2.114e+02 
 1.203e+00  4.594e+00  4.292e-01  1.755e+01  1.632e-01  1.640e+00 
 1.449e+01  5.532e+01  4.594e+00  2.114e+02  1.640e+00  1.755e+01 

My question is if I can give some "massage" to A before calling the LU 
decomposition? What I mean is something like a preconditioner which are 
often used when talking about iterative solvers. But when checking the 
documentation I couldn't find certain functionalities for the Matrix 
classes.
For me it would be great to know if the A matrix could be "rescued" somehow 
by certain solution techniques, preconditioning or whatever. Otherwise I 
have to think about the way of constructing A itself.

Thanks for giving me feedback!

Best
Simon

-- 
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/80e9dfe5-29ca-4414-9db2-5420caf3ab31n%40googlegroups.com.

[deal.II] triangulation.signals.post_refinement :connection of a member function

[Simon]

Dear all,

I am trying to modify my program in a way such that certain quantities are 
only updated if the underlying triangulation object changes during runtime 
(due to adaptivity). I want to call a specific member function of my Main 
Class 'Solid' if this happens. (Solid is the class which calls the typical 
functions like make_grid(), system_setup(),...)
The Triangulation class offers a interface to the 
boost::signals2 library. I guess that is predestinated for my purpose. 

Here is a snippet of my code:

//constructor of 'Solid'
:
dof_handler(triangulation),
/*many_more*/
qf_cell(degree+1),
n_q_points(qf_cell.size())
{
pointer_M = boost::make_shared>(/*constructor of 
 AnotherClass 
*/);
boost::signals2::connection tria_signal = 
triangulation.signals.post_refinement.connect( 
boost::signals2::signal::slot_type(&
*AnotherClass::mark_for_update*,  pointer_M .get()).track( pointer_M ) 
);
}

So the constructor of Solid initializes "pointer_M" which points to an 
object of "AnotherClass". The variable "tria_signal" in the line below is 
the one which causes the call of the member function &
*AnotherClass::mark_for_update. *
The above code compiles and actually calls *mark_for_update*  from 
*AnotherClass* . 

But what I'd like to achieve is that the signal triggers a member function 
of 'Solid', i.e. something like &*Solid::call_this_function. *So a 
naive approach is to simply replace &*AnotherClass::mark_for_update*  
by &*Solid::call_this_function* . (Both functions take no parameters 
and return void).
Unfortunately this change produces two error messages. 
(i) error: pointer to member type ‘void (Solid<2>::)()’ incompatible with 
object type ‘AnotherClass<2>’
(ii) error: return-statement with a value, in function returning 'void' 
[-fpermissive]

My understanding of the error message is that only member functions of 
'AnotherClass' can be connected to the signal. 
Is this true? 
And can I modify the parameters of the connect function somehow to actually 
connect a member function of Solid?

As always, any input is appreciated! :-)

Best
Simon

   

-- 
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/05d5ff07-e70a-412a-963e-83ab99efb11an%40googlegroups.com.

[deal.II] dealii in spack - approach to upgrade to latest version 9.3.3 via local git clone of spack

[Simon]

Dear all,

I am currently working with dealii version 9.1.1 and would like to install 
the latest version because I need access to the new element formulations. 
I cloned the spack repository and installed dealii via 
'spack install dealii'
and did all necessary instructions imposed by the dealii wiki entry on 
spack. However, I could not figure out from there how to install a newer 
version, for instance 9.3.3.

In my local working copy of the spack repository, I only have one master 
branch so far. I guess I first to need to upgrade spack somehow, in order 
to be able to install newer versions of a package, since 
'spack info dealii' 
only shows me versions up to 9.1.1, i.e. the one which I currently have 
installed.

1. So, should I first create a new branch in my local spack repo and then 
type 
'git fetch origin'
in order to bring spack up to date? 
I would like to keep version 9.1.1 as well, so I am a bit afraid of 
overwriting existing packages if I do everything on the master branch. Is 
this a valid approach at all or am I on the wrong track?
2. Is special care necessary, if a second version of dealii will be 
installed? 
As far as I know, spack is able to work with several versions of a package. 
However, I do not know how a second version of dealii affects all it's 
dependencies, for instance. UMFPACK. Will there newer versions installed 
'on-the-fly' as well, for the case they exist,...?

My lack of experience with spack makes me a bit careful to simple try 
things out, so I would appreciate to get some guidance how to accomplish 
this. 

Best
Simon






-- 
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/c1a94a60-9b7f-42c3-be36-64865b69791fn%40googlegroups.com.

[deal.II] finite element with shape functions having delta_ij property at quadrature points

[Simon]

Dear all,

I am solving a problem in 2d using FE_Q(2) elements and a gauss quadrature 
rule with (fe.degree +1) quadrature points in each co-ordinate direction, 
that is, I have in total nine quadrature points. My question pertains to 
the following:
At each cell, I need to approximate a field whose sampling (support) points 
are the quadrature points belonging to reduced integration, i.e, there are 
four quadrature points in my case and the four (shape) functions 
approximating my field should be designed as follows:
N_j (xi_k) = delta_{jk} ,
where xi_k are the coordinates of the quadrature points. So I need four 
(shape) functions, each of which is one at one of the four quadrature 
points and zero at the three others.  
That said, my ansatz is given by (the coefficents a(xi) are of course known)
f(xi) = a(xi_1) N_1(xi) + a(xi_2) N_2(xi) + a(xi_3) N_3(xi) + a(xi_4) N(xi)
and I need to evaluate the function f(xi) at the *nine* quadrature points.

What is the best way to do set up the FEValues object?
I have seen that there is a constructor for the FE_Q element which takes a 
Quadratute<1> object. I guess this would help me to define the (shape) 
functions pertaining to the field f(xi), but I think I can not evaluate 
this field at the nine quadrature points, because (i) their local 
coordinates are of course different in the new FEValues object and (ii) 
second, I would have to insert negative local coordinates for a set of 
them. 
Maybe I do not even need a FEValues object for my purpuse. As I said, I 
only need to do the approximation f(xi) and evaluate it at the nine 
quadrature points.

Any input is greatly appreciated!

Best
Simon



-- 
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/852af19c-dd9c-4c94-809f-e7b8b6df49d9n%40googlegroups.com.

[deal.II] cite specific content of dealii documentation

[Simon]

Dear all,

I wrote a thesis where my source code is based on dealii.

As described here,
https://www.dealii.org/publications.html,
I inserted the given bibtex entry in my document.
However, I also use specific content of
the documentation of the project function inside 
the vector tools namespace,
https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#ac6b404bf03cb2a742b290421cc2789fe,
in particular, the notes on the quadrature formula for the "mass matrix".

How do I cite this correctly? Simply listing all authors of the bibtex 
entry from above + the url?

Best
Simon

-- 
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/c6d23852-e7c2-4d27-a93a-75d6621225b7n%40googlegroups.com.

[deal.II] discretization of a analytical function on a second triangulation

[Simon]

Dear all:


I am solving a nonlinear PDE on a Triangulation T_1.

I know the analytical representation of a scalar function of two variables 
(no space co-ordinates, but two invariants of a quantity) 
and my goal is to find a discretized version of the analytical function (on 
a second, two-dimensional triangulation T_2).
The nodal values of T_2 should be the exact values corresponding to the 
analytical solution.
As for the values between the nodes, I want to start with a bilinear 
interpolation.

The coupling to my original PDE works as follows:
For each quadrature point defined on T_1, I compute the two input variables 
of the above function based on my solution vector, find the corresponding 
cell of T_2 and want to retrieve the interpolated nodal values. 

My idea is to use the function
VectorTools::interpolate ,
but I do not know a suitable Function object to hand over. 

Is this an appropriate approach at all?


Thank you,
Simon

-- 
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/14b1d8f6-f902-4414-a5c6-348ae3f5003en%40googlegroups.com.

[deal.II] shape function gradient at arbitrary points in the element

[Simon]

Dear all:

I have to compute the gradient with respect to real space co-ordinates of 
the shape functions belonging to a FE_Q element at arbitrary points in the 
element, that is, not only at quadrature points. 

What is the way to do this in dealii?
FEValues provides shape function values, gradients,... only at quadrature 
points. 

Thank you,
Simon

-- 
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/16afc8fa-bcfe-4a81-9ade-040c18f85596n%40googlegroups.com.

[deal.II] fourth-order referential deviatoric tensor Dev_P: why is 'S' used in the definition?

[Simon]

Dear all,


the fourth-order referential deviatoric tensor as returned by 
Physics::Elasticity::StandardTensors 
<https://www.dealii.org/current/doxygen/deal.II/classPhysics_1_1Elasticity_1_1StandardTensors.html><
 
dim >::Dev_P 
includes the fourth-order referential/spatial unit *symmetric* tensor 
Physics::Elasticity::StandardTensors 
<https://www.dealii.org/current/doxygen/deal.II/classPhysics_1_1Elasticity_1_1StandardTensors.html><
 
dim >::S = identity_tensor 
<https://www.dealii.org/current/doxygen/deal.II/symmetric__tensor_8h.html#ab3e890348aa219805e84f7d367e098c3>
().

In the literature, for instance G. A. Holzapfel: "Nonlinear solid 
mechanics. A Continuum Approach for Engineering" (2007), 
however, the general fourth-order unit tensor 
I_{ijkl} = delta_{ik} delta_{jl} is used to compute Dev_P. 

I implemented a hyperelastic material model with a volumetric / isochoric 
split of the strain energy function and it only converges when using 'S' - 
as dealii does it. Using the  general fourth-order unit tensor to define 
Dev_P,  my solver does not converge at all.  

Is it neccessary to use 'S' due to the way dealii stores and accesses the 
elements of symmetric tensors? 


Best
Simon


-- 
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/cc73ad68-f056-421d-b717-1bfcf1be9d4cn%40googlegroups.com.

[deal.II] interface between dealii and matlab for data transfer (using mex files?)

[Simon]

Dear all,

my question is specific and regards the data transfer between Matlab and a 
dealii program. In a nutshell, I have to call my dealii program from Matlab 
with a Matlab vector as argument (my dealii program needs the values in the 
vector). I figured out that Mex files are one way to transfer data to a c++ 
program.


Here is a snippet of my Matlab script:

function res = fun(p,y)
s = call_my_dealii_program_passing_the_vector_p
res = s - y;
end

My issue is *how to pass the vector 'p' to the dealii program*.

There are MEX files which provide an interface between Matlab and C++ 
programs. A standard MEX file looks like this:

//file mymex.cpp
#include "mex.h"
//include some more header files from c++ or dealii
//write here the whole class definition and implemenation

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
//this is the interface function to Matlab where also 'p' is passed
//in here, I will also create the instance of my dealii program to run the 
typical member //functions ike make_grid(), system_setup(),...
}

The above mex file is compiled within the Matlab gui using a c++ compiler 
(mex mymex.cpp). 

If I do not include any dealii headers in  mymex.cpp , the compilation 
works. 
However, including dealii headers results in a bunch of error messages when 
I compile the program in Matlab (mex mymex.cpp), because environment 
variables are not known to Matlab,...
Usually, I just call 'spack load dealii' followed by 'make run' to run my 
dealii programs, but I do not know how to 'forward' all this information to 
the Matlab compilation process.

All that said, has someone already worked with Mex files together with 
dealii?

I also appreciate other approaches to manage the data transfer between 
Matlab and dealii. 


Best
Simon

-- 
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/b047d35d-b3e3-465b-88d8-c708343086bfn%40googlegroups.com.

[deal.II] get the partition of the system matrix A associated with the unconstrained dofs

[Simon]

Dear all,


some degrees of freedom of my solution vector are constrained. 
When solsing the linear system A*x=b, I use the first approach as described 
in "Constraints on degrees of freedom". 
That is, I call
constraints.distribute_local_to_global(...,false)
and, after solving the linear system,
constraints.distribute(...). 

My question pertains to the following:
Say, I have in total 20 dofs and the dof with global dof index Zero is 
constrained. 
As a consequence, the first component of the rhs b is set to Zero as well as
the first row and column of the system matrix A (except the diagonal 
value). 
In a postprecessing step, I have to solve another linear system, however, 
with the system matrix being only the 19x19 matrix associated with the 19 
unconstrained dofs; I do not need the rhs anymore.

Is there a way to get this portion of the system matrix based on the 
existing system matrix (sparsity pattern)?


Thank you,
Simon

-- 
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/219e5a2e-06f3-400e-a3db-6e022020587an%40googlegroups.com.

[deal.II] visualizing a finite element field: different piecewise constant values within cells

[Simon]

Dear all,


my triangulation is two-dimensional and consists of quadrilaterals, i.e. 
four vertices per cell. Each cell has four quadrature points (qp's).

I visualize the data at the quadrature points (in paraview) by producing a 
discontinuous finite element field which has as many DoFs as my 
triangulation has qp's, that is, the qp values are somehow assigned to the 
nodes. 
In paraview (or other visualization programs) I would like to see four 
piecewice constant "colors" on every element: Sloppy speaking, each cell is 
divided into four subcells, and each of the subcells shows only one color 
which corresponds to the nodal value belonging to that subcell.

To this end, I think I need a finite element with piecewise constant shape 
functions. Each shape function is '1' at one quarter of the reference 
domain and '0' at the remaining three quarters. 
Or is there an easier way to do this?
If not, I probably would have to write my own finite element class, given 
that such an element does not exist in dealii?


Thank you,
Simon

-- 
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/0206f2a1-85cb-4970-9da8-d78de0bf6cc7n%40googlegroups.com.

[deal.II] automatic differentiation of pde solution wrt design variables in dealii

[Simon]

Dear all, 


I want to differentiate my pde solution U with respect to some design 
variables E, that is, computing the Jacobian dU/dE. 
I already implemented the Jacobian analytically, however, I would like to 
double check it using AD. 

" (In theory an entire program could be made differentiable. This could be 
useful in, for example, the sensitivity analysis of solutions with respect 
to input parameters. However, to date this has not been tested.)"

I found this information in the AD module. 
Also the relevant tuturials -- step33, step 70, step 71 -- do not cover the 
topic sensitivity analysis, but only AD at cell- and quadrature point 
level. 

That said, is it possible (with reasonable efforts) to make the entire 
dealii program AD differentiable, or is it recommendable to use other open 
source tools for that purpose?
My biggest concern is the solution of the linear system because the solver 
classes require in most cases a Vector and, consequently, the 
dependency u(E) can not be encoded.


Best
Simon

-- 
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/f9f35931-c2b1-4e49-9910-e142b7684804n%40googlegroups.com.

[deal.II] large condition number after solving linear systems

[Simon]

Dear all,

I am optimising parameters which are part of the pde I am working on.

To this end, I am solving n linear systems, where n is the number of 
parameters:
K * sensitivity[i] = rhs[i]i=1,...,n 
K is the (symmetric) tangent matrix which I already use to solve for the 
solution of my pde, that is, 
I just assemble the rhs vectors anew. 
The result is the Jacobian matrix 
J = [ sensitivity[1] ;...; sensitivity[n] ) 
which has as many rows as my pde has degrees of freedom. 
J has dimensions 320x8, that is, 
my pde has 320 dofs and I want to optimise 8 parameters. 
The solution of my pde (the dofs) are in the interval [0 to 1] and the 
parameters are in the interval [0 to 1e5] , i.e., the entries of J are 
rather small per nature (1e-6 to 1e-10). 
I evaluated the condition number of J in matlab: cond(J) = 1e13, which is 
too large given that the optimization algorithm works on J^t*J. 

Currently, I solve the above n linear systems using the direct solver 
UMFPACK. I know that there is no such thing like a preconditioner for 
direct solvers.
Thus, I solved the linear systems using the iterative SolverCG with a 
preconditioner:
const int solver_its = 1000
const double tol_sol = 1e-15
SolverControl solver_control(solver_its, tol_sol);
GrowingVectorMemory > GVM;
SolverCG > solver_CG(solver_control, GVM);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(tangent_matrix, 1.2);
//solve the linear systems...
However, the condition number of J is still around 1e13. 

My idea was to push the condition number of J at least some order of 
magnitudes when applying a preconditioner to my tangent matrix. 
Is this consideration reasonable?

Best
Simon

-- 
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/96c541fb-c5d2-4f90-8d9b-2a1f603d3430n%40googlegroups.com.

[deal.II] measuring cpu and wall time for assembly routine

[Simon]

Dear all,

I implemented two different versions to compute a stress for a given strain 
and want to compare the associated computation times in release mode.

version 1: stress = fun1(strain)  cpu time:  4.52  s  wall time:   
4.53 s
version 2: stress = fun2(strain) cpu time: 32.5s  wall time: 
167.5 s 

fun1 and fun2, respectively, are invoked for all quadrature points 
(1,286,144 in the above example) defined on the triangulation. My program 
is not parallelized.
In fun2, I call  find_active_cell_around_point 
<https://www.dealii.org/current/doxygen/deal.II/namespaceGridTools.html#a2e10aeb1c8e76110a84b6945eac3aaf0>
 
twice for two different points on two different (helper) triangulations and 
initialize two FEValues objects 
with the points ' ref_point_vol' and 'ref_point_dev' 
as returned by find_active_cell_around_point 
<https://www.dealii.org/current/doxygen/deal.II/namespaceGridTools.html#a2e10aeb1c8e76110a84b6945eac3aaf0>
 
.
FEValues<1> fe_vol(dof_handler_vol.get_fe(), 
Quadrature<1>(ref_point_vol),
update_gradients | update_values); 
   
FEValues<1> fe_values_energy_dev(this->dof_handler_dev.get_fe(), 
Quadrature<1>(ref_point_dev),
update_gradients | update_values); 
  

I figured out that the initialization of the two FEValues objects is the 
biggest portion of the above mentioned times.  In particular, if I comment 
the initialization out, I have 
cpu time: 6.54 s wall time: 6.55 s .

The triangulations associated with dof_handler_vol and dof_handler_dev are 
both 1d and store only 4 and 16 elements, respectively. That said, I am 
wondering why the initialization takes so long (roughly 100 seconds wall 
time in total) and why this causes a gap between the cpu and wall time.
Unfortunately, I have to reinitialize them anew whenever fun2 is called, 
because  the point 'ref_point_vol' (see Quadrature<1>(ref_point_vol)) is 
different in each call to fun2. 

Best
Simon



-- 
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/8f0e7843-c382-44bc-8e24-bbb65818c40cn%40googlegroups.com.

[deal.II] collect2 linking error when executing make

[Simon]

Dear all,

First of all, I wish you a happy and prosperous new year 2023!

I transferred my dealii project to a different computer and 
encounter an linking error (command make) after setting up 
the project (command cmake .).
The error message is attached.

What is the problem?

Best
Simon

-- 
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/ca0ce1f2-cba0-42e1-8ee4-2178806540ebn%40googlegroups.com.

[deal.II] possible to shift a support point (no vertex) of a triangulation?

[Simon]

Dear all,

I am approximating a scalar function F of one single argument by 
discretizing the argument over the interval [lb;ub]. 
Then, I interpolate the nodal values utilizing VectorTools::interpolate.
To generate the triangulation (a line), I call 
GridGenerator::hyber_cube(triangulation, lb, ub) and refine globally up to 
a desired level.
Either FE_Q(1) or FE_Q(2) elements are used. 

The issue is that there is one 'point' p in [lb;ub] for which the function 
must be Zero. 
In case of FE_Q(1) elements, I simply shifted the vertex which is closest 
to p, to p.
For FE_Q(2) elements, the closest point to p is a vertex *or* the support 
point in the center of the element. 
I think it is not possible to 'shift' a support point as the triangulation 
is made up of vertices only, right?

Best
Simon

-- 
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/e4c1dc91-7e1c-4317-93ca-ef8985f1313dn%40googlegroups.com.

[deal.II] B-spline basis functions to interpolate nodal values (control points) in 1d

[Simon]

Dear all,

my objective is to do an 1d-approximation 
F(x) = N_i(x) P_i 
where N_i are cubic B-spline basis functions and P_i known nodal values 
(control points). 
To assemble my linear system, I have to compute the first and second 
derivatives of F at the quadrature points of the triangulation. 
Since the number of control points P_i is quite moderate (less than 20), I 
think it is reasonable to evaluate all basis functions N_i in lieu of 
implementing an entire finite element class.

I scanned the dealii manual but could only find a CSpline class 
<https://www.dealii.org/current/doxygen/deal.II/classFunctions_1_1CSpline.html>
.
Just to make sure: 
B-Spline (basis) functions are not implemented in dealii yet, right?

So I probably will have to rely on an external library. 
Can someone make a library recommendation?
I read about tinyspline <https://github.com/msteinbeck/tinyspline> and 
splinelib <https://github.com/mfcats/SplineLib>, but have no knowledge 
regarding their performance, or if they are appropriate for usage in a 
dealii-program.

As I said, I mainly want to 
(i) interpolate control points with B-spline basis functions and
(ii) evaluate first and second derivatives at arbitrary points in the knot 
interval.

Best
Simon 

-- 
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/a26b043d-da5b-4ce6-80d5-c4d46734beefn%40googlegroups.com.

[deal.II] reference for citing the concept of vector-valued test functions?

[Simon]

Dear all,

is the deal.ii implementation of vector-valued test functions as described 
here 
  
also described in an article, paper, ... or is this a "deal.ii thing"?

I would like to make a reference in a publication where the concept of 
vector-valued test functions is explained in more detail, but do not want 
to cite a webpage.

Thank you!

-- 
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/5d15ce9c-1324-41bc-b832-2c578a31f093n%40googlegroups.com.

[deal.II] writing data to a file changes program output although this file is never used

[Simon]

Dear all,

my deal.ii program is called roughly 300 times from a matlab script for the 
purpose of a parameter identification. In essence, given some parameters, 
the matlab script calls deal.ii to compute a displacement field. I exchange 
the data transfer between deal.ii and matlab by reading from and writing 
into files. In particular, I store the parameters in a file with 13 digits 
after the decimal point since there are finite-difference calculations 
involved.
I double-checked that my simulations are deterministic, i.e. the parameters 
at all 300 calls to deal.ii are identical if I re-run the matlab script 
anew. By identical, I mean that all 13 digits coincide.

My deal.ii program has, besides the standard data structures 
(Triangulation, DoFHandler,...), a pointer to a class as member variable. 
This class has a member function, which is called during assembly of the 
linear system to calculate a stress tensor. Within this member function, I 
created a
std::ofstream out_data("/a/path/file.txt")
object and wrote the calculated stress into this file:
out_data

[deal.II] Re: Linking error when compiling examples: error: cannot find -lhdf5-shared

[Simon]

Is there already a solution to this issue?

I have the same problem as Amit Sharm;
hdf5 is loaded.

Best
Simon

amitshar...@gmail.com schrieb am Freitag, 31. März 2023 um 15:05:26 UTC+2:

> Dear Bruno,
>
> Yes, I loaded hdf5(version-1.10.7) and I tried same with latest version 
> available on spack.
>
> Thanks,
> Amit
>
> On Friday, March 31, 2023 at 6:27:13 PM UTC+5:30 bruno.t...@gmail.com 
> wrote:
>
>> Amit,
>>
>> Did you load the hdf5 package? You can check using `spack find --loaded`.
>>
>> Best,
>>
>> Bruno
>>
>> On Friday, March 31, 2023 at 4:51:02 AM UTC-4 amitshar...@gmail.com 
>> wrote:
>>
>>> Hello everyone,
>>>
>>> I installed dealii using spack manager. I started learning dealii using 
>>> examples but I am getting a linking error.
>>>
>>> Scanning dependencies of target step-1
>>> [ 33%] Building CXX object CMakeFiles/step-1.dir/step-1.cc.o
>>> [ 66%] Linking CXX executable step-1
>>> /home/amit/spack/opt/spack/linux-ubuntu20.04-skylake/gcc-9.4.0/binutils-2.38-e2nnae3w7g55p34q7m3tsngmmrol2iba/bin/ld.gold:
>>>  
>>> error: cannot find -lhdf5-shared
>>> collect2: error: ld returned 1 exit status
>>> make[3]: *** [CMakeFiles/step-1.dir/build.make:223: step-1] Error 1
>>> make[2]: *** [CMakeFiles/Makefile2:272: CMakeFiles/step-1.dir/all] Error 
>>> 2
>>> make[1]: *** [CMakeFiles/Makefile2:252: CMakeFiles/run.dir/rule] Error 2
>>> make: *** [Makefile:196: run] Error 2
>>>
>>> Thanks,
>>> Amit Sharma
>>>
>>>

-- 
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/74ca1a35-bc78-441c-8f44-408b61646fdcn%40googlegroups.com.

[deal.II] try catch of "ExcDereferenceInvalidCell"

[Simon]

Dear all,

I want to catch exceptions thrown by fePointEval.reinit(...), 
where fePointEval is an object of type FEPointEvaluation.

Here is a snippet of code demonstrating the issue:

try
{
 fePointEval.reinit(...);
}
catch(dealii::ExceptionBase & exception)
{
   std::cout<<"exception is catched"

[deal.II] inverse operation to AffineConstraints.condense() (recovering full system matrix)

[Simon]

Dear all, 

the "Constraints on Degrees of Freedom" module introduces several 
approaches to deal with constraints.
I implemented the "second" approach, that is, distributing the constraints 
on the fly when assembling the linear system using the 
constraints.distribute_local_to_global() function.
Since I plan to move from serial to distributed versions of my program in 
the future, 
it makes sense to stick to this approach as mentioned in the module.

After assembling and solving the linear system, I need the uncondensed 
system matrix and right hand side to perform a sensitivity analysis.
In particular, I need the rows and entries of the system matrix and RHS 
corresponding to constrained dofs. 
However, the constraints.distribute_local_to_global() function sets these 
entries all to zero (except the diagonal value of the system matrix).
My solution to deal with this is to basically copy my assembly function 
and, 
instead of distributing the constraints with the constraints object --
constraints.distribute_local_to_global() --
I distributed them with the cell iterator --
cell-> distribute_local_to_global(),
which does not care about the constraints.

This solution is of course very inefficient as for efficiency and code 
re-use.
As I said, it is basically a copy of my assembly function,
although most of the system matrix "is already there". 

So is there a more efficient way to recover the full system matrix based on 
the condensed system matrix?

Thanks for your help!

Best,
Simon


-- 
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/bcd7ffe7-34e0-4205-8296-c538a74b49efn%40googlegroups.com.

[deal.II] "is protected within this context" error for BlockSparseMatrix iterators

[Simon]

Dear all,

dealii::BlockSparseMatrix systemMatrix;
auto systemMatrixIteratorEnd = systemMatrix.end(0);
auto systemMatrixIteratorBegin = systemMatrix.begin(0);
if (systemMatrixIteratorBegin != systemMatrixIteratorEnd) std::cout<<
"Example"<::operator==(const 
dealii::BlockMatrixIterators::Accessor&) const 
[with BlockMatrixType = 
dealii::BlockMatrixBase >]’ is protected 
within this context
  184 |   return (accessor == other.accessor);

include/deal.II/lac/block_matrix_base.h:1475:3: note: declared protected 
here
 1475 |   Accessor::operator==(const Accessor &a) 
const


Everything works if BlockSparseMatrix is replaced by SparseMatrix.

Is there anything I miss here?

Best
Simon


-- 
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/b21d1915-c54a-4711-8b18-a8fa09d06a80n%40googlegroups.com.

[deal.II] anisotropic refinement with regard to global coordinate direction

[Simon]

Dear all,

attached is a minimal example in which I create two simple triangulations 
(tria1, tria2),
merge them, and finally anisotropically refine all cells 
(RefinementCase<3>::cut_xz).
I also write the merged triangulation to a file for visualizations sake.

My objective is to refine only in the global xy plane, that is, in global 
z-direction, I want to have one element everywhere (if hanging nodes allow 
that, of course).
For the given example,
RefinementCase<3>::cut_xz 
is correct for tria1, but not for tria2. 
I read the documenation which states that anisotropic refinement occurs 
w.r.t the local coordinate system of the cell, so I am not really surprised 
about the above result.

Anyway, is there a way to perform anisotropic with regard to the global 
system?
My ideas was to, based on the local orientation of the cell system, pass a 
different RefinementCase when marking the cells for refinement.
The TriaAccessor class offers the combined_face_orientation() function, but 
I am
not sure if it helps here.

Thank you for your input!

Best,
Simon


-- 
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/68df0485-a2dc-4c75-abdf-27dd52e87bb5n%40googlegroups.com.

#include 
#include 
#include 
#include 

#include 

using namespace dealii;

int main()
{

  Triangulation<3> tria, tria1, tria2;
  
  GridGenerator::hyper_cube_with_cylindrical_hole(tria1,
2.5,
7.5,
15.0,
1,
false);
  GridGenerator::hyper_cube(tria2, 7.5, 22.5);
  GridTools::shift(dealii::Point<3>{0,-15.0,-7.5}, tria2);
  // merge them
  dealii::GridGenerator::merge_triangulations<3>(tria1, tria2, tria);
  // execute global refinements
  for (unsigned int i = 0; i<2; ++i)
{
  for(const auto &cell : tria.active_cell_iterators())
cell->set_refine_flag(RefinementCase<3>::cut_xz);
  tria.execute_coarsening_and_refinement();
} 

  // output the mesh
  std::ofstream file("tria.vtk");
  dealii::GridOut   gridOut;
  dealii::GridOutFlags::Vtk vtkFlags(true, true, true, false);
  gridOut.set_flags(vtkFlags);
  gridOut.write_vtk(tria, file);

  return 0;
}

[deal.II] extract_surface_mesh() does not extract refined cells

[Simon]

Dear all,


attached is a minimal example which demonstrates that the function 
"extract_surface_mesh()" in the GridGenerator namespace 
does not extract refined cells. 
That is, it seems to work only for globally refined meshes. 
There is a comment in the implementation of this function:
// Create boundary mesh and mapping
// from only level(0) cells of the volume mesh.

Since the documentation states nothing about the usage of this function for 
adaptively
refined meshes (like it is done in flatten_triangulation, for instance), 
I am not sure whether this is desired behavior or not.

In my case, I probably do not need the full power of this function. 
All I want to do is to "copy" the solution values on a surface of the 
volume mesh 
to a surface mesh, and to associate a DoFHandler object to the surface mesh
that allows me to do further postprocessing.
So I do not care about the boundary and manifold id as well
as the return type, however, my volume mesh may contain hanging nodes.



Best regards,
Simon

-- 
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/f3252992-4761-41b1-a691-05c5265f89bfn%40googlegroups.com.

#include 
#include 
#include 
#include 

#include 

#include 

using namespace dealii;

int main()
{

  Triangulation<3> tria, triaTmp;
  
  dealii::Point<3> center;
  double radius = 2.5;
  unsigned int boundaryID = 2;

  GridGenerator::hyper_cube_with_cylindrical_hole(triaTmp,
radius,
7.5,
2.0,
1,
false);
  triaTmp.refine_global(2);

  // symmetry in global Z
  std::set::active_cell_iterator> cellsToRemove;
  for (const auto &cell : triaTmp.active_cell_iterators())
if ( cell->center()[2] > 1.0)
  cellsToRemove.insert(cell);
  dealii::GridGenerator::create_triangulation_with_removed_cells(triaTmp,
  cellsToRemove,
  tria);
  // one local refinement along hole
  for (const auto &cell : tria.active_cell_iterators())
for (unsigned int f = 0; f < dealii::GeometryInfo<3>::faces_per_cell; ++f)
  {
dealii::Point<3> faceCenter(cell->face(f)->center());
faceCenter[2] = 0.0;
if ((std::abs(center.distance(faceCenter)) < radius))
  cell->set_refine_flag();
  }
  tria.execute_coarsening_and_refinement();

  // set boundary ID for surface mesh
  for (const auto &cell : tria.active_cell_iterators())
for (unsigned int f = 0; f < dealii::GeometryInfo<3>::faces_per_cell; ++f)
  if(std::abs(cell->face(f)->center()[2]) < 1e-5 && cell->face(f)->at_boundary())
cell->face(f)->set_boundary_id(boundaryID);

  // extract surface mesh
  dealii::Triangulation<2, 3> tria2d;
  dealii::GridGenerator::extract_boundary_mesh(
tria,
tria2d,
std::set{boundaryID});

  // output the volume mesh
  std::ofstream file("tria.vtk");
  dealii::GridOut   gridOut;
  dealii::GridOutFlags::Vtk vtkFlags(true, true, true, false);
  gridOut.set_flags(vtkFlags);
  gridOut.write_vtk(tria, file);

  // and the surface mesh
  std::ofstream file2d("tria2d.vtk");
  gridOut.write_vtk(tria2d, file2d);

  return 0;
}

[deal.II] construct TrilinosWrappers::SparseMatrix from TrilinosWrappers::BlockSparseMatrix to use Amesos direct solver

[Simon]

Dear all,

there are several discussions, e.g.
https://groups.google.com/g/dealii/c/KU7g_YCNVRE/m/h7ab-2gRCAAJ?hl=en 
or 
https://groups.google.com/g/dealii/c/DdMiTjuPx6Q/m/UeH4fmmBAwAJ?hl=en
about the usage of TrilinosWrappers::BlockSparseMatrix in the context of 
direct solvers 
in a distributed framework.
Still, only TrilinosWrappers::SparseMatrix can be feeded to 
TrilinosWrappers::SparseDirect solvers.  

I am solving a 2x2 block system with iterative and direct solvers (decision 
is made at run-time via .prm file). 
Clearly, I do not need a BlockSparseMatrix for a direct solve, however, 
to have a similar code base also for the iterative solve (where I need a 
block preconditioner), 
I decided to use block quantities for each solver type. 

It is not a big deal to copy a BlockVector into a Vector, but what is an 
efficient scheme (my problem is non-linear) to 
copy a BlockSparseMatrix into a SparseMatrix as needed for the direct 
solvers coming with trilinos?
Conceptually, I am looking for a constructor 
SparseMatrix (const BlockSparseMatrix &) 
or a copy_from function. 
For instance, the std::memcpy function is used in 
SparseMatrix::reinit (const Epetra_CrsMatrix & *, **copy_values* = true) ,
but I am not sure whether it can be used when translating between Block and 
standard SparseMatrix. 

Can you think of an elegant solution or must the non-zero values be copied 
via a loop
over the individual blocks?

Best regards,
Simon



-- 
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/70bcce5b-73c4-4550-ac0c-368ce549e0bfn%40googlegroups.com.

Re: [deal.II] construct TrilinosWrappers::SparseMatrix from TrilinosWrappers::BlockSparseMatrix to use Amesos direct solver

[Simon]

 "This. I don't see a different way of doing things. Though I might just 
always 
create both objects and then call distribute_local_to_global() with one or 
the 
other matrix object, depending on your run-time flags. "

Your approach sounds more reasonable than copying data back and forth at 
every Newton step
for a possibly large linear system.

I use a BlockSolver for the solution and right hand side of the linear 
system to compute the various norms (non-linear solver) 
of the the individual blocks.
For a direct solver, I would now use a SparseMatrix as matrix type and a 
BlockVector as vector type.
At each solve, I convert the BlockVector to a standard Vector, which I 
guess is not too expensive.
This mix (SparseMatrix and BlockVector) should not cause any undue 
problems, right?
Or would you not recommend this?

Best
Simon





Wolfgang Bangerth schrieb am Montag, 18. Dezember 2023 um 18:03:13 UTC+1:

> On 12/18/23 04:47, Simon wrote:
> > 
> > Can you think of an elegant solution or must the non-zero values be 
> copied via 
> > a loop
> > over the individual blocks?
>
> This. I don't see a different way of doing things. Though I might just 
> always 
> create both objects and then call distribute_local_to_global() with one or 
> the 
> other matrix object, depending on your run-time flags.
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/9b61d0e7-6ea3-4db9-aa6a-578304bee319n%40googlegroups.com.

Re: [deal.II] construct TrilinosWrappers::SparseMatrix from TrilinosWrappers::BlockSparseMatrix to use Amesos direct solver

[Simon]

Timo,

" What I have done in several codes (and what is used in ASPECT), is to 
always create block matrices and vectors but to put all DoFs into a single 
block when a direct solver is used. This decision can be made at rintime. 
Inside the solve routine you can just access block_matrix.block(0,0), which 
is a SparseMatrix. Same for the vector. This way, there are no different 
code paths, different types, or copies involved. "

In terms of readability, this sounds like the best solution. 

My motivation to use block structures even for direct solvers is to compute
norms of the pressure or the velocity, if we talk about the Stokes system.
Or in general, having access to the individual blocks to write them to 
files,...

If I am not mistaken, you create a nested finite element system with one 
block

FESystem<3>feVelocity( FE_Q<3>(degree), 3 );
FE_DGQ<3>   fePressure(degree-1);
FESystem<3>nestedFE ( feVelocity, 1, fePressure, 1 ); 

resulting in a 1x1 block system.
But it is not clear to me how to compute the norms of velocity and pressure 
in that case.
So far, I used a Block indices object as returned by 
DoFTools::count_dofs_per_fe_block()
since I compute the l2-norm only of the unconstrained DoFs. 

That said, how do you access the velocity and pressure block when working 
with direct solvers?

Thank you,
Simon
Timo Heister schrieb am Dienstag, 19. Dezember 2023 um 20:02:30 UTC+1:

> What I have done in several codes (and what is used in ASPECT), is to 
> always create block matrices and vectors but to put all DoFs into a single 
> block when a direct solver is used. This decision can be made at rintime. 
> Inside the solve routine you can just access block_matrix.block(0,0), which 
> is a SparseMatrix. Same for the vector. This way, there are no different 
> code paths, different types, or copies involved.
>
> --
> Timo Heister
> http://www.math.clemson.edu/~heister/
> --
> *From:* dea...@googlegroups.com  on behalf of 
> Wolfgang Bangerth 
> *Sent:* Monday, December 18, 2023 4:37:34 PM
> *To:* dea...@googlegroups.com 
> *Subject:* Re: [deal.II] construct TrilinosWrappers::SparseMatrix from 
> TrilinosWrappers::BlockSparseMatrix to use Amesos direct solver 
>  
> This Message Is From An External Sender: Use caution when opening links or 
> attachments if you do not recognize the sender.
>
>
> On 12/18/23 12:44, Simon wrote:
> >
> > I use a BlockSolver for the solution and right hand side of the
> > linear system to compute the various norms (non-linear solver) of the
> > the individual blocks. For a direct solver, I would now use a
> > SparseMatrix as matrix type and a BlockVector as vector type. At each
> > solve, I convert the BlockVector to a standard Vector, which I guess
> > is not too expensive. This mix (SparseMatrix and BlockVector) should
> > not cause any undue problems, right? Or would you not recommend
> > this?
> No. This seems like a reasonable approach to me.
> Best
>   W.
>
> --
> The deal.II project is located at 
> https://nam12.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=05%7C02%7Cheister%40clemson.edu%7Cb26c1543ea3f4a770d1a08dc00118fcd%7C0c9bf8f6ccad4b87818d49026938aa97%7C0%7C0%7C638385322953674220%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000%7C%7C%7C&sdata=FCQhyskq54Y%2B9Uc7oyauIeXMsI5CLzO12RE%2FCDG1w3Y%3D&reserved=0
>  
> <http://www.dealii.org/>
> For mailing list/forum options, see 
> https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=05%7C02%7Cheister%40clemson.edu%7Cb26c1543ea3f4a770d1a08dc00118fcd%7C0c9bf8f6ccad4b87818d49026938aa97%7C0%7C0%7C638385322953674220%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000%7C%7C%7C&sdata=QKshKRkLcezx3WBxZW4FajigvudlynOE6GClmS9nNnY%3D&reserved=0
>  
> <https://groups.google.com/d/forum/dealii?hl=en>
> ---
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit 
> https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2F15994792-761a-4899-b7e6-af2fc6ba2932%2540colostate.edu&data=05%7C02%7Cheister%40clemson.edu%7Cb26c1543ea3f4a770d1a08dc00118fcd%7C0c9bf8f6ccad4b87818d49026938aa97%7C0%7C0%7C638385322953674220%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000%7C%7C%7C&sdata=Lmjub869SyEbNQ5DTAN4RNnzG41%2B2XEWpPs6s97c9OU%3D&reserve

[deal.II] access all elements of a row of a trilinos sparse matrix

[Simon]

Dear all,


I want to build the scalar product between rows of a
TrilinosWrappers::SparseMatrix and a TrilinosWrappers::MPI::Vector. 
I implemented this by looping over the elements of the rows via iterators 
(it = systemMatrix.begin(row) , itend = systemMatrix.end(row) 
-> while (it != itend) { ... it++; } )
Inside while, I calculate it->value (to get the matrix element) times 
vec(it->column()) 
(to get the corresponding vector element) and sum it up.
I calculate the scalar product for a set of rows and store the result in a 
FullMatrix 
which I finally write to a file from the master process,
after gathering the local full matrices from all processors.
Attached is a small example which demonstrates this in practice.

In a parallel program, vec(it->column()) only works if it->column() is a 
locally owned element of the vector. So I wrapped this read access in 
if( vec.in_local_range(it->column() ) .
Another problem is the following:
Say row 5 lives on processor 0 and has non-zero entries at columns 
10,15,20,25,30.
Columns 10,15,20 live on processor 0 as well but columns 25,30 live on 
processor 1. 
Then, I will never add the contributions to the scalar product from columns 
25 and 30 because on processor 1, 
systemMatrix.begin(5) == systemMatrix.end(5)
which is also stated in the docu.

How can I deal with this?


Best,
Simon

-- 
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/03085575-2931-4755-abcd-7b7a31df102dn%40googlegroups.com.

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

#include 
#include 
#include 

#include 

#include 
#include 
#include 

#include 

#include 
#include 
#include 
#include 

namespace dealiiLA
{
  using namespace dealii::TrilinosWrappers;
} // namespace dealiiLA

int main(int argc, char *argv[])
{

const unsigned intmasterProcessRank = 0;
dealii::Utilities::MPI::MPI_InitFinalize mpiInitialization(argc, argv, 1);
MPI_Comm const &mpiCommunicator(MPI_COMM_WORLD);
unsigned int processorRank = 
dealii::Utilities::MPI::this_mpi_process(mpiCommunicator);
dealii::ConditionalOStream pCout(std::cout, processorRank == 
masterProcessRank);

// make grid
dealii::parallel::distributed::Triangulation<3> tria(mpiCommunicator);
dealii::GridGenerator::hyper_cube(tria);
tria.refine_global(2);

// distribute dofs
dealii::DoFHandler<3> dofHandler(tria);
dofHandler.distribute_dofs(dealii::FE_Q<3>(1));

dealii::IndexSet locallyOwnedDofs = dofHandler.locally_owned_dofs();
dealii::IndexSet locallyRelevantDofs;
dealii::DoFTools::extract_locally_relevant_dofs(dofHandler, 
locallyRelevantDofs);

dealii::AffineConstraints constraintMatrix;
constraintMatrix.clear();
constraintMatrix.reinit(locallyRelevantDofs);
dealii::DoFTools::make_hanging_node_constraints(dofHandler, 
constraintMatrix);
constraintMatrix.close();

// sparse matrix and sparsity pattern
dealiiLA::SparsityPattern sparsityPattern;
sparsityPattern.reinit(locallyOwnedDofs, mpiCommunicator);
dealii::DoFTools::make_sparsity_pattern(dofHandler,
  sparsityPattern,
  constraintMatrix,
  false,
  
dealii::Utilities::MPI::this_mpi_process(
mpiCommunicator));
sparsityPattern.compress();

dealiiLA::SparseMatrix systemMatrix;
systemMatrix.reinit(sparsityPattern);

// create dummy vector
dealii::TrilinosWrappers::MPI::Vector vec;
vec.reinit(locallyOwnedDofs, mpiCommunicator);
for(unsigned int i=0; i rows{10, 40, 70, 100, 120}; //nDoFs = 125
dealii::FullMatrix matrix(rows.size(), 1);
// populate fullmatrix
for(unsigned int i = 0; i < rows.size(); ++i)
{
  unsigned int idx = rows[i];
  auto it = systemMatrix.begin(idx);
  auto itend = systemMatrix.end(idx);
  // loop over row idx (first contribution)
  while(it != itend)
  {
matrix(i, 0) += vec.in_local_range((*it).column()) ? ((*it).value() + 
1) * vec((*it).column()) : 0.0;
it++;
  }
  // second contribution
  if(vec.in_local_range(idx))
matrix(i, 0) += vec(idx);
}

// send local matrices to processor 0
std::vector> gathered = 
  dealii::Utilities::MPI::gather(mpiCommunicator, matrix, 0);

// write the matrix from processor 0
if(processorRan

[deal.II] trilinos sparsity pattern: end iterator of last locally owned row triggers trilinos error

[Simon]

Dear all,

attached is a minimal example with 125 dofs, where dofs with indices 0-74 
live on proccesor 0 and dofs with indices 75-124 live on processor 1
(two MPI processes).

Consider the following code to loop over the rows of a Trilinos sparsity 
pattern:
// make_sparsity_pattern,...
for(unsigned int row = 0; row < dofHandler.n_dofs(); ++row)
{
   if( ! sparsityPattern.row_is_stored_locally(row)) continue;
   auto it = sparsityPattern.begin(row); 
   auto itend = sparsityPattern.end(row); 
   // do something with the columns
   while(it != itend) { it->column();...}
}

For row=74 on processor 0, there is an error 

dealii::TrilinosWrappers::SparsityPatternIterators::Accessor::visit_present_row()
"An error with error number -1 occurred while calling a Trilinos function"

The mentioned Trilinos function is ExtractGlobalRowCopy() and the error 
number
suggests an access to a row which is not locally owned. 
This makes sense to me because 
sparsityPattern.end(74) calls visit_present_row(75),
which is indeed not locally owned.

Not sure whether I miss something or this is a bug.
My expectation was to use the begin and end iterators on all locally owned 
rows.

Thank you,
Simon

-- 
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/cb270b3a-8c45-40cf-9efe-2a579fdab28an%40googlegroups.com.

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

#include 
#include 

#include 

#include 
#include 
#include 

#include 

#include 

namespace dealiiLA
{
  using namespace dealii::TrilinosWrappers;
} 

int main(int argc, char *argv[])
{
dealii::Utilities::MPI::MPI_InitFinalize mpiInitialization(argc, argv, 1);
MPI_Comm const &mpiCommunicator(MPI_COMM_WORLD);

dealii::parallel::distributed::Triangulation<3> tria(mpiCommunicator);
dealii::GridGenerator::hyper_cube(tria);
tria.refine_global(2);

dealii::DoFHandler<3> dofHandler(tria);
dealii::FE_Q<3> fe(1);
dofHandler.distribute_dofs(fe);

dealii::DoFRenumbering::component_wise(dofHandler);

dealii::IndexSet locallyOwnedDofs = dofHandler.locally_owned_dofs();

dealiiLA::SparsityPattern sparsityPattern;
sparsityPattern.reinit(locallyOwnedDofs, mpiCommunicator);
dealii::DoFTools::make_sparsity_pattern(dofHandler,
  sparsityPattern,
  dealii::AffineConstraints(),
  false,
  
dealii::Utilities::MPI::this_mpi_process(
mpiCommunicator));
sparsityPattern.compress();

for(unsigned int i=0; i

Re: [deal.II] trilinos sparsity pattern: end iterator of last locally owned row triggers trilinos error

[Simon]

"Would you be willing to create a complete example, and open a bug report 
for
this on github?"

Sure. I opened an issue 
https://github.com/dealii/dealii/issues/16414

Do you think this bug has a quick solution?
I am interested in the columns of the rows. 
I did not see a different way to do this than using the row iterators.

Best,
Simon

On Thursday, January 4, 2024 at 9:50:36 PM UTC+1 Wolfgang Bangerth wrote:


Simon: 

> attached is a minimal example with 125 dofs, where dofs with indices 0-74 
> live on proccesor 0 and dofs with indices 75-124 live on processor 1 
> (two MPI processes). 
> 
> Consider the following code to loop over the rows of a Trilinos sparsity 
pattern: 
> // make_sparsity_pattern,... 
> for(unsigned int row = 0; row < dofHandler.n_dofs(); ++row) 
> { 
>if( ! sparsityPattern.row_is_stored_locally(row)) continue; 
>auto it = sparsityPattern.begin(row); 
>auto itend = sparsityPattern.end(row); 
>// do something with the columns 
>while(it != itend) { it->column();...} 
> } 
> 
> For row=74 on processor 0, there is an error 
> 
> 
dealii::TrilinosWrappers::SparsityPatternIterators::Accessor::visit_present_row()
 

> "An error with error number -1 occurred while calling a Trilinos 
function" 
> 
> The mentioned Trilinos function is ExtractGlobalRowCopy() and the error 
number 
> suggests an access to a row which is not locally owned. 
> This makes sense to me because 
> sparsityPattern.end(74) calls visit_present_row(75), 
> which is indeed not locally owned. 
> 
> Not sure whether I miss something or this is a bug. 

A bug. And a timely one in view of the fact that I have been thinking about 
exactly this issue recently: 
https://github.com/dealii/dealii/pull/16406 
Would you be willing to create a complete example, and open a bug report 
for 
this on github? 

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4f17ab8a-c292-4601-83dc-16e43a9c69a9n%40googlegroups.com.

[deal.II] mapping of global dof indices after solution transfer

[Simon]

Dear all,

I am storing a std::vector of size N holding global DOF 
indices living at some nodes of my mesh where we have measured experimental 
values.
After doing a SolutionTransfer motivated by AMR, I want to have the new 
global DOF indices at the same nodes, so the new vector of indices has the 
same size N.

I feel the only way to do this is to store the indices along with their 
support points before I enter the AMR cycle (so far I only store the 
indices) and determine the new indices after refinement by looping over the 
new mesh and check the support points for equality.

Is there a more efficient solution for this problem?

Thanks,
Simon

-- 
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/0e4febd6-d884-4584-a33c-49ba8f255c37n%40googlegroups.com.

[deal.II] "constraints for degree of freedom ... should not be stored ... " after SolutionTransfer

[Simon]

Dear all,

I am doing a SolutionTransfer in the same fashion as described in the docu 
of SolutionTransfer class. 
After calling
-triangulation.execute_coarsening_and_refinement(); 
I proceed with
-dofHandler.distribute_dofs(feSystem);
-constraints.clear();
-DoFTools::make_hanging_node_constraints(dofHandler, constraints);
where my feSystem has 2 blocks (3 copies of FE_Q<3(1) and 1 copy of 
FE_DGQ<3>(0))

The function 
DoFTools::make_hanging_node_constraints(dofHandler, constraints)

throws

  dealii::types::global_dof_index dealii::AffineConstraints<
number>::calculate_line_index(size_type) const [with number = double; 
dealii::types::global_dof_index = unsigned int; size_type = unsigned int]
The violated condition was: 
local_lines.is_element(line_n)
Additional information: 
The index set given to this constraints object indicates constraints
for degree of freedom 167 should not be stored by this object, but a
constraint is being added. 

Before refinement, dofHandler stored 138 DoFs, after refinement 202. 

I could not reproduce the error on a minimal example.
However, it is not clear to me why
make_hanging_constraints() reports an error associated with dof index 167, 
if I cleared the constraints object before hand and can be certain that my 
DoFHandler has 202 DoFs. 
So the sequence of distributing DoFs, followed by clearing constraints and 
calling make hanging node constraints should work independently of how my 
program arrived at these lines of code. Is this not true?

Any general hints as for the problem source and narrowing down the problem? 

Thanks,
Simon
 
 



-- 
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/8c7a9d55-549e-4750-9141-88add9f91470n%40googlegroups.com.

Re: [deal.II] "constraints for degree of freedom ... should not be stored ... " after SolutionTransfer

[Simon]

No, the program that produced the error is serial not parallel.

Best,
Simon

Wolfgang Bangerth schrieb am Montag, 10. Juni 2024 um 04:30:33 UTC+2:

> On 6/8/24 07:51, Simon wrote:
> > 
> > Before refinement, dofHandler stored 138 DoFs, after refinement 202.
> > 
> > I could not reproduce the error on a minimal example.
> > However, it is not clear to me why
> > make_hanging_constraints() reports an error associated with dof index 
> 167,
> > if I cleared the constraints object before hand and can be certain that 
> my 
> > DoFHandler has 202 DoFs.
> > So the sequence of distributing DoFs, followed by clearing constraints 
> and 
> > calling make hanging node constraints should work independently of how 
> my 
> > program arrived at these lines of code. Is this not true?
> > 
> > Any general hints as for the problem source and narrowing down the 
> problem?
>
> Did you initialize the AffineConstraints object with index sets for the 
> locally owned and/or locally relevant DoFs? If so, you will have to 
> re-initialize it with the now larger/different number of DoFs you have 
> after 
> mesh refinement.
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/66d6a3fd-eb61-45ab-9ca1-509ed97bdc39n%40googlegroups.com.

Re: [deal.II] "constraints for degree of freedom ... should not be stored ... " after SolutionTransfer

[Simon]

I figured out the problem. 

In my make_constraints() function, I re-initialized the AffineConstraints 
object indeed with 
the locally releveant DoFs as you said (although this is not necessary in 
the serial version of my code). 

Thanks for pointing out to this!

-Simon

Wolfgang Bangerth schrieb am Montag, 10. Juni 2024 um 19:05:01 UTC+2:

>
> On 6/10/24 00:13, Simon wrote:
> > No, the program that produced the error is serial not parallel.
>
> I think then I'd need to see the whole program :-)
>
> Best
> W.
>

-- 
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/11ee8fef-a9e5-45e2-968c-e8382b743c3bn%40googlegroups.com.

[deal.II] Handling NaN Values in TrilinosWrappers::SparseMatrix with Amesos_SuperluDist

[Simon]

Dear all,

In the attached example, I deliberately assign NaN values to a 
TrilinosWrappers::SparseMatrix to test how the Amesos_SuperluDist direct 
solver handles such inputs. The code is run using deal.II version 9.4.0.

Unfortunately, the program results in a segmentation fault. Specifically, 
no exception (TrilinosWrappers::SolverDirect::ExcTrilinosError) is thrown 
when the matrix containing NaN values is passed to the solver's solve 
function.

In contrast, when running a sequential program with SparseDirectUMFPACK, 
the solver correctly handles the issue, and I am able to catch the 
exception SparseDirectUMFPACK::ExcUMFPACKError.

Can this be considered a Trilinos bug? Has this issue been addressed in 
newer versions of deal.II or Trilinos?

Clearly, assigning NaN values to a matrix does not make much sense. 
However, in the optimized release version of my program, bad inputs can 
occasionally lead to cases where the cell matrix contains NaN values. In 
such situations, I rely on catching exceptions (e.g., from 
SparseDirectUMFPACK) to handle the error gracefully.

Best,
Simon

-- 
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 visit 
https://groups.google.com/d/msgid/dealii/8a68096d-c0ba-4edc-976c-17788993c18an%40googlegroups.com.
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 

#include 
#include 

#include 
#include 

#include 
#include 

int main(int argc, char *argv[])
{
// Initialize MPI
dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

const MPI_Comm mpi_communicator = MPI_COMM_WORLD;

// Create a 1D triangulation with 2 cells
dealii::Triangulation<1> triangulation;
dealii::GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(1); // Create 2 cells

// Finite element and DoF handler
const dealii::FE_Q<1> fe(1); // Linear elements
dealii::DoFHandler<1> dof_handler(triangulation);
dof_handler.distribute_dofs(fe);

// IndexSet for locally owned DoFs
dealii::IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
dealii::IndexSet locally_relevant_dofs;
dealii::DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);

// Trilinos sparse matrix
dealii::TrilinosWrappers::SparseMatrix system_matrix;
dealii::DynamicSparsityPattern dsp(locally_owned_dofs);
dealii::DoFTools::make_sparsity_pattern(dof_handler, dsp, dealii::AffineConstraints{}, true);
system_matrix.reinit(locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator);
dealii::SparsityTools::distribute_sparsity_pattern(dsp,
 locally_owned_dofs,
 mpi_communicator,
 locally_relevant_dofs);

// Local matrix
const unsigned int dofs_per_cell = fe.dofs_per_cell;
dealii::FullMatrix local_matrix(dofs_per_cell, dofs_per_cell);
std::vector local_dof_indices(dofs_per_cell);

// Assign NaN to local matrix and distribute it to global system matrix
for (const auto &cell : dof_handler.active_cell_iterators())
{
  cell->get_dof_indices(local_dof_indices);
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
  local_matrix(i, j) = std::numeric_limits::quiet_NaN(); // Fill local matrix with NaN values
  cell->distribute_local_to_global(local_matrix, system_matrix);
}

// Finalize the global matrix assembly
system_matrix.compress(dealii::VectorOperation::add);

// Right-hand-side vector and solution vector
dealii::TrilinosWrappers::MPI::Vector rhs(locally_owned_dofs, mpi_communicator);
dealii::TrilinosWrappers::MPI::Vector solution(locally_owned_dofs, mpi_communicator);

// Solver
dealii::TrilinosWrappers::SolverDirect::AdditionalData additionalData;
additionalData.solver_type = "Amesos_Superludist";

dealii::SolverControl  solverControl;
dealii::TrilinosWrappers::SolverDirect directSolver(solverControl, additionalData);

// Attempt to solve the system
directSolver.solve(system_matrix, solution, rhs);

// Output the solution
solution.print(std::cout);

  return 0;
}

Re: [deal.II] Handling NaN Values in TrilinosWrappers::SparseMatrix with Amesos_SuperluDist

[Simon]

*If you run the program in a debugger, can you find out where the segfault 
happens?*

Yes, the starting point of the trace is 
"TrilinosWrappers::SolverDirect::do_solve", and the segmentation fault 
occurs in the libsuperlu_dist library. 

On reviewing the do_solve function, it seems the exception 
TrilinosWrappers::SolverDirect::ExcTrilinosError is expected to be thrown. 
While Amesos_Mumps solver correctly throws this exception, unfortunately, 
Amesos_superludist does not. :(

As you mentioned, it might be worthwhile to try a newer version of deal.II, 
which provides wrappers for the more advanced Amesos2 and Tpetra packages.

*Perhaps a better approach than seeing whether a linear solver succeeds 
would
be to first check whether the matrix/rhs have NaNs before you hand them off 
to
the solver. *

My motivation for avoiding such checks was to keep the optimized release 
code free of
unnecessary overhead, trusting that the solver would throw appropriate 
exceptions when feeded
with bad input. However, I probably can not rely on that for all solvers, 
and manually checking for NaN
is the more reliable approach.

-Simon
On Wednesday, December 18, 2024 at 11:53:42 PM UTC+1 Wolfgang Bangerth 
wrote:


Simon: 

> Unfortunately, the program results in a segmentation fault. Specifically, 
no 
> exception (TrilinosWrappers::SolverDirect::ExcTrilinosError) is thrown 
when 
> the matrix containing NaN values is passed to the solver's solve 
function. 

If you run the program in a debugger, can you find out where the segfault 
happens? 


> In contrast, when running a sequential program with SparseDirectUMFPACK, 
the 
> solver correctly handles the issue, and I am able to catch the exception 
> SparseDirectUMFPACK::ExcUMFPACKError. 
> 
> Can this be considered a Trilinos bug? Has this issue been addressed in 
newer 
> versions of deal.II or Trilinos? 

Segfaults should never happen. Better error messages are always better, but 
I 
don't know whether this has been addressed. Perhaps this is your chance to 
try 
a newer version of deal.II? ;-) 


> Clearly, assigning NaN values to a matrix does not make much sense. 
However, 
> in the optimized release version of my program, bad inputs can 
occasionally 
> lead to cases where the cell matrix contains NaN values. In such 
situations, I 
> rely on catching exceptions (e.g., from SparseDirectUMFPACK) to handle 
the 
> error gracefully. 

Perhaps a better approach than seeing whether a linear solver succeeds 
would 
be to first check whether the matrix/rhs have NaNs before you hand them off 
to 
the solver. This can be done cheaply via the Vector::l1_norm() and 
SparseMatrix::frobenius_norm() functions that simply add up the (squares) 
or 
entries. If the result is NaN, you know something is amiss. 

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.
To view this discussion visit 
https://groups.google.com/d/msgid/dealii/cecde726-6a28-4a14-95c9-395c19250a26n%40googlegroups.com.

[deal.II] Re: Decoupling FECollection and QCollection

[Simon Sticko]

Hi,
When you call reinit on hp::FEValues you can send in which of the elements 
in the collection you want to reinitialize with.
That is, when you call reinit you call it as 

fe_values.reinit(cell, q_index, mapping_index, fe_index);

where the additional integers specify which element in the collections you 
want to use.

/Simon

>

-- 
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: Decoupling FECollection and QCollection

[Simon Sticko]

Well, I guess it's a question of how one should interpret this sentence. 
Anyway, if you want to avoid calling reinit with extra parameters during 
assembling (might be simpler anyway) you can probably achieve the same 
thing by adding the same element twice to FECollection.
For  example:

fe_collection.push_back(FE_Q(fe_order1));
fe_collection.push_back(FE_Q(fe_order2));
fe_collection.push_back(FE_Q(fe_order2));

q_collection.push_back(QGauss(quad_order1);
q_collection.push_back(QGauss(quad_order2);
q_collection.push_back(QGauss(quad_order3);

Best
/Simon

-- 
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: fe_enriched and step-47

[Simon Sticko]

Hi.
Regarding this:

> On Thursday, April 6, 2017 at 1:51:41 PM UTC+2, Denis Davydov wrote:
> 
> The issue with integrating and quadrature would still remain.


I've been implementing an algorithm to generate bulk and surface quadrature 
given a level set function as boundary description. It's mainly an 
implementation of this algorithm: dx.doi.org/10.1137/140966290. 
I'm hoping to finish this within the coming weeks and contribute with this 
to dealii.

Best 
Simon

-- 
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] 2d meshes and orientation

[Konrad Simon]

Dear deal.ii community,

I stumbled over something in 2d meshes that confuses me a bit. Supoose I 
have 2 squares (a left and a right one) sharing one edge. If everything is 
fine all lines and faces have standard orientation, i.e., 
line_orientation==true and face_orientation==true for all lines and faces. 
Note that lines and faces are NOT the same geometric entity here since 
their orientation is determined differently (faces have normals and lines 
have tangents).

Now, if I rotate the right square by 90 degrees I get (on the shared edge) 
non-matching face normals but all line orientations do match.

If I rotate the right square by 180 degrees I get (on the shared edge) 
non-matching face normals and non-matching line orientations.

If I rotate the right square by 270 degrees I get (on the shared edge) 
matching face normals and non-matching line orientations.

This is at least my intuition. When I query these information for each 
rotation case I see that my intuition is confirmed for lines but the face 
orientations are always true.

Is this a bug?

I am asking because it seems that 2d RaviartThomas elements (which make use 
of face normals) are broken for some edge cases (3d RT I could fix).

Best,
Konrad

-- 
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/f739c450-c37d-4d24-bf5b-2778dec5d501n%40googlegroups.com.

Re: [deal.II] 2d meshes and orientation

[Konrad Simon]

Dear Wolfgang,

On Tuesday, April 6, 2021 at 5:53:15 AM UTC+2 Wolfgang Bangerth wrote:

> It's possible this is indeed a bug. In most cases, we run the 2d meshes 
> through the mesh orienter, so we rarely see 2d meshes that are not 
> correctly 
> oriented and it wouldn't surprise me if we have the kind of bug you 
> describe 
> where things almost but not quite work correctly.


I see. Many thanks for the input. This is my impression as well. These 
things rarely happen in 2d.
One example would be simple periodic meshes which cause line orientations 
not to match (this case is caught). In more
complicated scenarios involving periodicity and rotations of faces it can 
happen though that lines do match but normals do not.


> The questions then are (i) where does this need to be fixed, and (ii) how 
> many 
> tests fail if we apply the fix? One of the things I learned (and I think 
> you 
> did as well) is that for bugs in these kinds of layers, it is quite 
> possible 
> that one ends up with 100 failing tests and then it takes a *lot* of work 
> to 
> look through all of these to figure out whether the new output is correct 
> or 
> not. But then, as mentioned above, it's also possible that only one or two 
> tests fail because we almost never see these kinds of meshes.
>

Sound worth trying. Please see also this new PR 
 which is work in progress.

-- 
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/0125f6dc-b61c-4db6-84da-b101b6759f1bn%40googlegroups.com.

[deal.II] Re: interpolate FE_Nelelec

[simon...@gmail.com]

Hi,
I don't have any experience with FE_Nedelec, but the error message states 
that the function you are trying to interpolate has 1 component, while the 
element has dim components.

If you implemented your own Function, you need to make sure it has 
dim-components by calling the constructor of Function with dim as 
inargument:

template 
class CustomFunction : public Function
{
public:
CustomFunction()
: Function(dim) // The function has dim components
{}

double
value(const Point & point,
const unsigned int component = 0) const override
{
// Return value will depend on component here
return 0;
}
};


Best,
Simon

On Friday, April 16, 2021 at 5:07:01 PM UTC+2 johnsmi...@gmail.com wrote:

> Hello,
>
> It seems I am unable to find a function similar to 
> VectorTools::interpolate for FE_Nedelec finite elements. The existing 
> implementation of this functions gives the following run-time error:
>
> *An error occurred in line <556> of file 
> 
> */numerics/vector_tools.templates.h> in function*
>
> * void dealii::VectorTools::interpolate(const dealii::Mapping spacedim>&,*
>
> * const DoFHandlerType&, const dealii::Function typename*
>
> * VectorType::value_type>&, VectorType&, const dealii::ComponentMask&) 
> [with int *
>
> *dim = 3; int spacedim = 3; VectorType = dealii::Vector; 
> DoFHandlerType =*
>
> * dealii::DoFHandler; typename VectorType::value_type = double]*
>
> *The violated condition was: *
>
> * dof_handler.get_fe().n_components() == function.n_components*
>
> *Additional information: *
>
> * Dimension 3 not equal to 1.*
>
>
> Looks like VectorTools::interpolate does not like the vector-valued finite 
> elements. Could you please point me in the right direction?
>
>
> Best, 
>
> John
>

-- 
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/6d051c9e-8840-4e5d-8021-f8007a1dc543n%40googlegroups.com.

Re: [deal.II] Re: interpolate FE_Nelelec

[Konrad Simon]

Dear John,

I had a look at the pdf you sent. I noticed some conceptual inconsistencies 
that are important for the discretization. Let me start like this: The 
curl-curl-problem that you are trying to solve is - according to your 
description of the gauge - actually a curl-curl + grad-div problem and 
hence a Laplace-problem that describes a Helmholtz-decomposition. In other 
words in order to get a well-posed simulation problem you need more 
boundary conditions. You already noticed that since your curl-curl-matrix 
is singular (the kernel is quite complicated and contains all longitudinal 
waves). 

You have of course the option to solve (0.0.2) with a Krylov-solver but 
then you need to make sure that the right-hand side is in the orthogonal 
complement of this kernel before each iteration which is quite difficult. I 
would not recommend that option.

Another point is that if you have a solution for your field A like in 
(0.0.4) you can not have a similar representation for T. This is 
mathematically not possible.

Since you not care about the gauge let me tell you how I would tackle this: 
The reason  is  - to come back to my original point - you are missing a 
boundary condition that makes you gauge unique. Since you apply natural 
boundary conditions (BCs) on A you must do so as well to determine div(A). 
This second BC for A is applied to a different part of the system that is 
usually neglected in the literature and A is partly determined from this 
part (this part then describes the missing transversal waves). The 
conditions you mention on curl(A) are some what contradictory to the space 
in which A lives. Nedelec elements which you must use (you can not use FE_Q 
to enforce the conditions) cannot generate your desired T_0.

There are some principles when discretizing these problems (which are not 
obvious) that you MUST adhere to (choice of finite elements, boundary 
conditions, what is the exact system etc) if you want to get a stable 
solution and these are only understood recently. I am solving very similar 
problems (with Deal.II) in a fluid context and will be happy to share my 
experiences with you. Just email me: konrad.si...@uni-hamburg.de

Best,
Konrad

On Monday, April 19, 2021 at 5:51:17 AM UTC+2 Wolfgang Bangerth wrote:

> On 4/18/21 12:24 PM, John Smith wrote:
> > 
> > Thank you for your reply. It is a bit difficult to read formulas in 
> text. So I 
> > have put a few questions I have in a pdf file. Formulas are better 
> there. It 
> > is attached to this message. May I ask you to have a look at it?
>
> John:
>
> If I understand your first question right, then you are given a vector 
> field J 
> and you are looking for a vector field T so that
> curl T = J
> I don't really have anything to offer to this kind of question. There are 
> many 
> vector fields T that satisfy this because of the large null space of curl. 
> You 
> have a number of condition you would like to "approximate" but there are 
> many 
> ways to "approximate" something. In essence, you have two goals: To 
> satisfy 
> the equation above and to approximate something. You have to come up with 
> a 
> way to weigh these two goals. For example, you could look to minimize
> F(T) = ||curl T - J||^2 + alpha ||T-T_desired||^2
> where T_desired is the right hand side of (0.0.5) and you have to 
> determine 
> what alpha should be.
>
> As for your other questions: (0.0.9) is indeed called "interpolation"
>
> The issue with the Nedelec element (as opposed to a Q(k)^dim) field is 
> that 
> for the Nedelec element,
> \vec phi(x_j)
> is not independent of
> \vec phi(x_k)
> and so you can't choose the matrix in (0.0.8) as the identity matrix. You 
> realize this in (0.0.13). The point I was trying to make is that (0.0.13) 
> cannot be exactly satisfied if T(x_j) is not a vector parallel to \hat 
> e_j, 
> which in general it will not be. You have to come up with a different 
> notion 
> of what it means to solve (0.0.13) because you cannot expect the left and 
> right hand side to be equal.
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/75104172-4893-4ae6-af8e-8ec2c12d335an%40googlegroups.com.

Re: [deal.II] interpolate FE_Nelelec

[Konrad Simon]

Dear John,

As Jean-Paul pointed out the entries of the solution vector for a Nedelec 
element have a different meaning. The entries are weights that result from 
edge moments (and in case of higher order also face moments and volume 
moments), i.e, integrals on edges. Nedelec elements have a deep geometric 
meaning: they discretize quantities that you can integrate (a physicist 
would say "measure") along 1D-lines in 3D.

Btw, if you are using Nedelec in 2d be aware that it has some difficulties 
on complicated meshes. In 3D the lowest order is fine. NedelecSZ (for 2D 
lowest order and 3D all orders), as Jean-Paul pointed out, are another 
option. They are meant to by-pass certain mesh orientation issues on 
complicated geometries (I can tell you a bit about that since currently I 
am chasing some problems there).

Best,
Konrad

On Tuesday, April 20, 2021 at 9:23:23 PM UTC+2 Jean-Paul Pelteret wrote:

> Hi John,
>
> Unfortunately, you’ve fallen into the trap of confusing what the entries 
> in the solution vector mean for the different element types. For Nedelec 
> elements, they really are coefficients of a polynomial function, so you 
> can’t simply set each coefficient to 1 to visualise the shape function 
> (like one might be inclined to try for Lagrange type polynomials). If you 
> want a little piece of code that will plot the Nedelec shape functions for 
> you, then take a look over here:
> https://github.com/dealii/dealii/issues/8634#issuecomment-595715677
>
> BTW. In case you’re not aware, a newer (and probably more robust) 
> implementation of a Nedelec element is the FE_NedelecSZ 
> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> that 
> uses hierarchical shape functions to form the higher order bases. The PhD 
> thesis of S. Zaglmayr that describes the element is mentioned in the class 
> documentation.  
>
> Best,
> Jean-Paul
>
> On 20. Apr 2021, at 15:47, John Smith  wrote:
>
>
> Dear Konrad,
>
> Thank you for your help! I have recently modified the step-3 tutorial 
> program to save FE_Nedelec<2> fe(0) shape function into *.csv files. It is 
> attached to this message. The result of this program is not what I would 
> expect it to be. I have also attached a pdf file which compares my 
> expectations (marked as “Expectations” in the pdf file) with the results of 
> the modified step-3 program (marked as “Reality” in the pdf file). My 
> expectations are based on the literature I have. My question is: why the 
> shape functions of FE_Nedelec<2> fe(0) differ from my expectations? I guess 
> my expectations are all wrong. I would also appreciate if you will share 
> with me a reference to a paper or a book, so I can extend my modest library.
>
> Thanks!
>
> John
>
> P.S. I have sent you an e-mail as you have requested. Could you please 
> check your spam folder? 
>
>
>
> On Monday, April 19, 2021 at 8:47:29 AM UTC+2 Konrad Simon wrote:
>
>> Dear John,
>>
>> I had a look at the pdf you sent. I noticed some conceptual 
>> inconsistencies that are important for the discretization. Let me start 
>> like this: The curl-curl-problem that you are trying to solve is - 
>> according to your description of the gauge - actually a curl-curl + 
>> grad-div problem and hence a Laplace-problem that describes a 
>> Helmholtz-decomposition. In other words in order to get a well-posed 
>> simulation problem you need more boundary conditions. You already noticed 
>> that since your curl-curl-matrix is singular (the kernel is quite 
>> complicated and contains all longitudinal waves). 
>>
>> You have of course the option to solve (0.0.2) with a Krylov-solver but 
>> then you need to make sure that the right-hand side is in the orthogonal 
>> complement of this kernel before each iteration which is quite difficult. I 
>> would not recommend that option.
>>
>> Another point is that if you have a solution for your field A like in 
>> (0.0.4) you can not have a similar representation for T. This is 
>> mathematically not possible.
>>
>> Since you not care about the gauge let me tell you how I would tackle 
>> this: The reason  is  - to come back to my original point - you are missing 
>> a boundary condition that makes you gauge unique. Since you apply natural 
>> boundary conditions (BCs) on A you must do so as well to determine div(A). 
>> This second BC for A is applied to a different part of the system that is 
>> usually neglected in the literature and A is partly determined from this 
>> part (this part then describes the missing transversal waves). The 
>> conditions you mention on curl(A) are some what contradictory to the space 
>> in which 

Re: [deal.II] interpolate FE_Nelelec

[Konrad Simon]

Hi John,

Maybe I can give some hints of how I would approach this problem (these are 
just some quick thoughts):
Write this problem as a vector Laplace problem

curl(nu*curl(A)) - grad(div(A)) = J

which you then have to write as a system of PDEs. Note, this can only be 
the same as

curl(nu*curl(A)) = J

if J is a curl itself.

If you do not want to impose conditions on the divergence then you have two 
options:

1.) You use sigma = div(A) as an auxiliary variable and impose natural 
boundary conditions (i.e., the ones that you get through integration by 
parts) on A*n and nu*curl(A)xn, where n is the outward unit normal. This 
means that sigma is a 0-form and A is a 1-form -> sigma is then discretized 
(for example) by FE_Q(n+1) and A by FE_Nedelec(n)

System:

(A1) sigma + div(A) = 0
(B1) grad(sigma) + curl(nu*curl(A)) = J

For the weak form you multiply (A1) with a FE_Q test function and integrate 
the second summand by parts. Secondly, you multiply (B1) with a Nedelec 
test function and also integrate the second summand by parts. You will see 
your natural boundary conditions pop up.

2.) You use sigma = nu*curl(A) as an auxiliary variable and impose 
essential boundary conditions on A*n and nu*curl(A)xn, where n is the 
outward unit normal. This means that sigma is interpreted as a 1-form and A 
is a 2-form -> sigma is then discretized (for example) by FE_Nedelec(n) and 
A by FE_RaviartThomas(n)

System:

(A2) (1/nu)*sigma - curl(A) = 0
(B2) curl(sigma) -  grad(div(A)) = J

For the weak form you multiply (A2) with a Nedelec test function and 
integrate the second summand by parts. Secondly, you multiply (B2)
with a Raviart-Thomas test function and also integrate the second summand 
by parts. You will NOT see your natural boundary conditions pop up since 
conditions on A*n and nu*curl(A)xn are essential in this case. You need to 
enforce them on the function spaces directly. In deal.ii you do so by using 
project_boundary_values_curl_conforming_l2 

 
or project_boundary_values_div_conforming 


Note that the additional boundary conditions make the system invertible and 
if J is a curl it will turn out that div(A)=0.

There is a field in mathematics that is called finite element exterior 
calculus (FEEC) that answers the question of stability when solving such 
problems. See this book 
. I guess Chapters 
4, 5 and 6 are most interesting for you.

> I do not understand when a geometry gets complicated. Is a toroid inside a 
> sphere, both centred at the origin, a complicated geometry? To start with, 
> I will use a relatively simple mesh. I will not use local refinement, only 
> global. I can do without refinement as well, i.e. making the mesh in gmsh.
>
> Yes, I would like to know more about orientations issues on complicated 
> geometries. Could please tell me about it? I would very much prefer to use 
> FE_Nedelec without hierarchical shape functions. The first order 3D 
> FE_Nedelec will do nicely provided the orientation issues will be 
> irrelevant to the mesh.
>
 Lowest order Nedelec and all other elements I mentioned are fine on the 
meshes you need I guess.

Hope that helps.
Best,
Konrad  

-- 
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/c2773841-4e03-42a1-9f1e-e1a41c0933c9n%40googlegroups.com.

[deal.II] Re: Accessing cell data from its ID

[simon...@gmail.com]

Hi,
I'm not sure I understand exactly what you want to do, but if you know 
which cell you want in a triangulation by level and index you can get its 
center by creating the cell iterator. For example:

const int dim = 2;
Triangulation triangulation;
GridGenerator::hyper_cube(triangulation);

const unsigned int cell_level = 0, cell_index = 0;
const auto cell = Triangulation::active_cell_iterator(&triangulation,

  
cell_level,

  
cell_index);
std::cout << cell->center() << std::endl;


Maybe this answer is useful too:
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#can-i-convert-triangulation-cell-iterators-to-dofhandler-cell-iterators

Best,
Simon

On Thursday, April 29, 2021 at 10:48:30 AM UTC+2 corbin@gmail.com wrote:

> Hi all,
>
> I'm looking to map DoF from the bottom of an extruded volumetric mesh to 
> the corresponding DOF on the surface of the mesh. Right now, I'm saving the 
> cell "number"/ID of each cell in the bottom layer of the volume mesh, and 
> want to match it with a surface cell by comparing their centers
>
>- I could loop over all the volume cells until I find the one matching 
>the surface cell center, but as I have the cell numbers, I would rather 
>access the cell center from the cell id (which looks to me like an 
> unsigned 
>int)
>- something like get_cell(cell_id)->center  --- is this possible?
>- Alternatively, I was thinking that a periodic boundary condition 
>linking the bottom of the volume mesh with the top might be able to 
> achieve 
>the same thing, but I wasn't sure how to go about setting it up
>
> Any pointers would be appreciated!
> Corbin
>

-- 
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/3833b90d-1855-494f-953c-38fa2f4ce868n%40googlegroups.com.

Re: [deal.II] GridTools::find_cells_adjacent_to_vertex: how to avoid multiple calling ?

[Simon Wiesheier]

Hello Wolfgang,

I actually have a Time Measurement in my code implemented.
For my current computations the difference (one call vs two calls per
vertex) is neglible. But my model will become much bigger and therefore I'd
like to avoid unnecessary calls to that function.

Since I can pass only one dof_handler to that function, I guess there is no
"obvious" way to make only one call, is it?

-Simon

Wolfgang Bangerth  schrieb am Fr., 30. Apr. 2021,
18:16:

> On 4/30/21 10:13 AM, Simon wrote:
> >
> > auto cells_ref =
> GridTools::find_cells_adjacent_to_vertex(dof_handler_ref,
> > /*vertex*/);
> > auto cells_tmp =
> GridTools::find_cells_adjacent_to_vertex(dof_handler_tmp,
> > /*vertex*/);
> >
> > I am not sure how "expensive" this function actually is. I have to call
> this
> > function for all vertices of my triangulation, so currently I call it
> > "2*n_vertices()" times.
> >
> > I´d like to call this function only once for a given vertex, i.e.
> > "n_vertices()"times.
> > Is there a way to realize this in combination with two different
> dof_handlers?
>
> There are probably ways to work around this, but I'd like to repeat a
> point
> Bruno made the other day on this very mailing list:
>
> "Unless you have measured that this is a bottleneck in you code, you
> should
> use what's the more readable. If there is a difference between these two
> codes, it would need to be in a hot loop to matter. My advice is to write
> easy
> to understand code and once you have made sure that the code works, then
> use a
> profiler to find the bottleneck and optimize the code."
>
> I think this is good advice!
>
> Best
>   W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/966f38ec-5ce2-f879-c9c4-2df21ab1110b%40colostate.edu
> .
>

-- 
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/CAM50jEtCWSTAhzEi7uPmgvkC94DMGJnP43AAhOuJw_97KCVT7g%40mail.gmail.com.

Re: [deal.II] GridTools::find_cells_adjacent_to_vertex: how to avoid multiple calling ?

[Simon Wiesheier]

Hi Jean-Paul,

your suggestion worked fine, thanks a lot!

Best
Simon

Am Sa., 1. Mai 2021 um 07:20 Uhr schrieb Jean-Paul Pelteret <
jppelte...@gmail.com>:

> Hi Simon,
>
> You could call GridTools:: find_cells_adjacent_to_vertex() with the
> triangulation instead of a DoFHandler, and then just convert the returned
> iterators to the type used by each DoFHandler using the method outlined
> here:
>
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#can-i-convert-triangulation-cell-iterators-to-dofhandler-cell-iterators
>
> Best,
> Jean-Paul
>
> On 30. Apr 2021, at 18:37, Simon Wiesheier 
> wrote:
>
> Hello Wolfgang,
>
> I actually have a Time Measurement in my code implemented.
> For my current computations the difference (one call vs two calls per
> vertex) is neglible. But my model will become much bigger and therefore I'd
> like to avoid unnecessary calls to that function.
>
> Since I can pass only one dof_handler to that function, I guess there is
> no "obvious" way to make only one call, is it?
>
> -Simon
>
> Wolfgang Bangerth  schrieb am Fr., 30. Apr. 2021,
> 18:16:
>
>> On 4/30/21 10:13 AM, Simon wrote:
>> >
>> > auto cells_ref =
>> GridTools::find_cells_adjacent_to_vertex(dof_handler_ref,
>> > /*vertex*/);
>> > auto cells_tmp =
>> GridTools::find_cells_adjacent_to_vertex(dof_handler_tmp,
>> > /*vertex*/);
>> >
>> > I am not sure how "expensive" this function actually is. I have to call
>> this
>> > function for all vertices of my triangulation, so currently I call it
>> > "2*n_vertices()" times.
>> >
>> > I´d like to call this function only once for a given vertex, i.e.
>> > "n_vertices()"times.
>> > Is there a way to realize this in combination with two different
>> dof_handlers?
>>
>> There are probably ways to work around this, but I'd like to repeat a
>> point
>> Bruno made the other day on this very mailing list:
>>
>> "Unless you have measured that this is a bottleneck in you code, you
>> should
>> use what's the more readable. If there is a difference between these two
>> codes, it would need to be in a hot loop to matter. My advice is to write
>> easy
>> to understand code and once you have made sure that the code works, then
>> use a
>> profiler to find the bottleneck and optimize the code."
>>
>> I think this is good advice!
>>
>> Best
>>   W.
>>
>> --
>> 
>> Wolfgang Bangerth  email: bange...@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.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/966f38ec-5ce2-f879-c9c4-2df21ab1110b%40colostate.edu
>> .
>>
>
> --
> 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/CAM50jEtCWSTAhzEi7uPmgvkC94DMGJnP43AAhOuJw_97KCVT7g%40mail.gmail.com
> <https://groups.google.com/d/msgid/dealii/CAM50jEtCWSTAhzEi7uPmgvkC94DMGJnP43AAhOuJw_97KCVT7g%40mail.gmail.com?utm_medium=email&utm_source=footer>
> .
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/9301657C-CE72-4129-A738-A7AAA52F0C0F%40gmail.com
> <https://groups.google.com/d/msgid/dealii/9301657C-CE72-4129-A738-A7AAA52F0C0F%40gmail.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEu90yajFfKoN9G0wCRTNOoQFAt73BckneoGTfDhpCKhsQ%40mail.gmail.com.

Re: [deal.II] GridTools::find_cells_adjacent_to_vertex: how to avoid multiple calling ?

[Simon Wiesheier]

Thanks for all answers!

I've got a similar question:

I looked up the Implementation of the mentioned function. Basically the
function loops over all active cells and compares those global vertices
with the one passed as second argument, so far so good.

But is there a similar function which accepts a cell_iterator (the given
vertex is located on that cell)  so that inside the function body only the
neigbor and the neighbor's neighbor are possible cells on which the given
vertex can be located?
So I'd like to tell the function where it should begin its "search", in
order to avoid looping over all cells (worst case).

There is a function which is doing this for a patch of neighbor_cells
(element-patch), but I didn't find a function for the problem description
from above.

Is such a function not yet implemented or is there already one available?

Best
Simon

luca.heltai  schrieb am Mo., 3. Mai 2021, 09:47:

> Dear Simon,
>
> one solution you have is to pass a *Triangulation*, and then generate the
> dof cell accessors using the constructor
>
> typename DoFHandler::cell_iterator cell1(*tria_iterator, dh1);
> typename DoFHandler::cell_iterator cell2(*tria_iterator, dh2);
>
>
> Best,
> Luca.
>
> > On 30 Apr 2021, at 18:37, Simon Wiesheier 
> wrote:
> >
> > Hello Wolfgang,
> >
> > I actually have a Time Measurement in my code implemented.
> > For my current computations the difference (one call vs two calls per
> vertex) is neglible. But my model will become much bigger and therefore I'd
> like to avoid unnecessary calls to that function.
> >
> > Since I can pass only one dof_handler to that function, I guess there is
> no "obvious" way to make only one call, is it?
> >
> > -Simon
> >
> > Wolfgang Bangerth  schrieb am Fr., 30. Apr.
> 2021, 18:16:
> > On 4/30/21 10:13 AM, Simon wrote:
> > >
> > > auto cells_ref =
> GridTools::find_cells_adjacent_to_vertex(dof_handler_ref,
> > > /*vertex*/);
> > > auto cells_tmp =
> GridTools::find_cells_adjacent_to_vertex(dof_handler_tmp,
> > > /*vertex*/);
> > >
> > > I am not sure how "expensive" this function actually is. I have to
> call this
> > > function for all vertices of my triangulation, so currently I call it
> > > "2*n_vertices()" times.
> > >
> > > I´d like to call this function only once for a given vertex, i.e.
> > > "n_vertices()"times.
> > > Is there a way to realize this in combination with two different
> dof_handlers?
> >
> > There are probably ways to work around this, but I'd like to repeat a
> point
> > Bruno made the other day on this very mailing list:
> >
> > "Unless you have measured that this is a bottleneck in you code, you
> should
> > use what's the more readable. If there is a difference between these two
> > codes, it would need to be in a hot loop to matter. My advice is to
> write easy
> > to understand code and once you have made sure that the code works, then
> use a
> > profiler to find the bottleneck and optimize the code."
> >
> > I think this is good advice!
> >
> > Best
> >   W.
> >
> > --
> > 
> > Wolfgang Bangerth  email: bange...@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.
> > To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/966f38ec-5ce2-f879-c9c4-2df21ab1110b%40colostate.edu
> .
> >
> > --
> > 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/CAM50jEtCWSTAhzEi7uPmgvkC94DMGJnP43AAhOuJw_97KCVT7g%40mail.gmail.com
> .
>
&

Re: [deal.II] GridTools::find_cells_adjacent_to_vertex: how to avoid multiple calling ?

[Simon Wiesheier]

Thanks to both of you!

The implementation of a own function was the best solution to my issue. In
there I basically ask the neighbor-cells and again those neighbors if the
given vertex is located on that cells. My function is a bit faster than
your suggestion Luca, but not that much.

Nonetheless, the GridTools::Cache-Class is actually much faster than
calling GridTools::adjacent_cell_to_vertex for all vertices of my
triangulation. So thanks for that hint, I guess I can you use some other
features of that class in my program.

Best
Simon

Am Mi., 5. Mai 2021 um 14:01 Uhr schrieb luca.heltai :

> In addition to that, I’d suggest you use the
> GridTools::Cache::get_vertex_to_cell_map class, which builds the map once,
> and stores it for future reference. This one uses
> GridTools::vertex_to_cell_map, which builds the map with one pass, instead
> of looping over all cells.
>
> Indeed, if you want to call this only for one vertex, your choice is
> better. But if you need this for several vertices, I’d build the map once,
> and then use the map all subsequent times.
>
>
> https://www.dealii.org/current/doxygen/deal.II/classGridTools_1_1Cache.html#a1e51292fbaeae9d7d3fcec2f9eb6a8dd
>
>
> Best,
> Luca.
>
> > On 4 May 2021, at 22:35, Wolfgang Bangerth 
> wrote:
> >
> > On 5/4/21 12:15 PM, Simon Wiesheier wrote:
> >> I looked up the Implementation of the mentioned function. Basically the
> function loops over all active cells and compares those global vertices
> with the one passed as second argument, so far so good.
> >> But is there a similar function which accepts a cell_iterator (the
> given vertex is located on that cell)  so that inside the function body
> only the neigbor and the neighbor's neighbor are possible cells on which
> the given vertex can be located?
> >> So I'd like to tell the function where it should begin its "search", in
> order to avoid looping over all cells (worst case).
> >> There is a function which is doing this for a patch of neighbor_cells
> (element-patch), but I didn't find a function for the problem description
> from above.
> >> Is such a function not yet implemented or is there already one
> available?
> >
> > I don't think so, but it shouldn't be very hard to add an argument to
> the function that by default is simply a default-constructed cell iterator
> indicating "no initial guess given" but if not at its default value
> indicates to start search on this cell and its neighbors.
> >
> > Best
> > W.
> >
> > --
> > 
> > Wolfgang Bangerth  email: bange...@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.
> > To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/47cadef0-3287-db4f-27e4-a4923757374a%40colostate.edu
> .
>
> --
> 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/EE1CEADB-D4D8-4E2F-8263-1F21416AF1AF%40gmail.com
> .
>

-- 
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/CAM50jEsrrW_ykF8m_wH63jH%2B0a7v-xHMDkJvnf6Wa-aqV0OC8A%40mail.gmail.com.

[deal.II] p4est user pointer and distributed::triangulation

[Konrad Simon]

Dear all, 

I am currently writing an interface to some p4est functions. The structure 
of p4est makes it often necessary to pass data around through a user 
pointer of whatever type (void*). When creating a new triangulation this 
pointer is set to „this“ (the triangulation itself), see for example

https://www.dealii.org/current/doxygen/deal.II/distributed_2tria_8cc_source.html#l2989

Is this pointer used in other interfaces somehow? p4est itself does not 
touch it so I am wondering if I can reset it to anything I‘d like.

Any Ideas?

Best,
Konrad

-- 
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/f23950cf-7688-45d2-b986-140ccc3a58f9n%40googlegroups.com.

Re: [deal.II] p4est user pointer and distributed::triangulation

[Konrad Simon]

Thank you, Wolfgang! This is exactly what I am looking for. I don‘t need to 
access the pointer from outside, I am trying to extend the p4est-interface 
itself.

Best,
Konrad

On Friday, May 14, 2021 at 12:42:48 AM UTC+2 Wolfgang Bangerth wrote:

> On 5/13/21 3:05 PM, Konrad Simon wrote:
> > 
> > I am currently writing an interface to some p4est functions. The 
> structure of 
> > p4est makes it often necessary to pass data around through a user 
> pointer of 
> > whatever type (void*). When creating a new triangulation this pointer is 
> set 
> > to „this“ (the triangulation itself), see for example
> > 
> > 
> https://www.dealii.org/current/doxygen/deal.II/distributed_2tria_8cc_source.html#l2989
>  
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fdistributed_2tria_8cc_source.html%23l2989&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cb4cadd28da84445fd29308d91652e2c3%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637565367510208994%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=v2G100jGfI6SbGwV35raSFjMxfnJMJurP2F2qAc61Ls%3D&reserved=0
> >
> > 
> > Is this pointer used in other interfaces somehow? p4est itself does not 
> touch 
> > it so I am wondering if I can reset it to anything I‘d like.
>
> Konrad,
> we use this user pointer in many of the functions that are called back 
> from 
> p4est. Take a look at
> RefineAndCoarsenList::refine_callback
> for example (line 761 of the file), and many other functions that you can 
> find 
> by searching for
> forest->user_pointer
> in this file.
>
> We consider the p4est forest object as internal to the p::d::T class, so I 
> would expect that you can't access it, or expect that you can use any of 
> its 
> content for anything other than what the p::d::T class does with it.
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/7fbd1adf-c6d4-4656-ab6e-7c833e5e92cen%40googlegroups.com.

Re: [deal.II] projectin of gradients of the discrete solution to nodes - two different discontinuous approaches

[Simon Wiesheier]

Thank you!

I failed to see that I used the FE_DGPMonomial instead of FE_DGQ.
With the latter it makes as you said acutally no difference.

However I am wondering why FE_DGPMonomial gives my quite good results
as well. As I said, for "fine meshes" the results of the local and global
approach are nearly the same and for coarse meshes, i.e. just one or two
times globally refined meshes, the biggest deviation also reaches not more
than 5%.
I didn´t have to use FE_DGPMonomial so far. I´ve seen that for degree=1
this element has only 3 dofs per cell for a scalar problem in 2D, i.e I can
not extrapolate the 4 qp-values to the nodes. But let´s assume I use this
element for my local approach.
-What is the relations of the resulting 3 nodel values to the 4 qp-values?
Is it still something like minimizing the squared difference or is the idea
a different one?
-And why is there a difference at all between the local and global
approach? Should the global approach not result in a block-diagonal matrix
as well?

I couldn´t figure out the purpose of the disontinuous legendre-polynomials
by having a look at the documentation.

-Simon

Am Mo., 24. Mai 2021 um 17:54 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:

> On 5/22/21 1:22 AM, Simon wrote:
> >
> > I am aware that I solve a linear system in the first case whereas in the
> > second case not, but my idea was the following: The first approach
> minimizes
> > the squared difference, but if I have as many dofs_per_cell as qps then
> of
> > course the squared difference can be minimized to zero for each qp, i.e.
> the
> > qp values can directly be transferred to the nodes. And if my
> understanding is
> > correct, this is exactly what the second approach does. So in the first
> > approach I do not do it "directly" but due to its definition and the
> number of
> > dofs it should do the same as second.
> > Of course I compared this approaches in my program. However depending on
> the
> > mesh size there are deviations up to percent, with finer meshes this
> > difference reduces.
> > Can this deviation be argumented away with the standard argument "this
> is the
> > numerics..." or is there is a mathematical difference and both
> approaches do
> > something different?
>
> For discontinuous elements, the mass matrix you compute is block diagonal
> and
> the projection can be computed cell-by-cell -- so it shouldn't make a
> difference whether you solve the global problem or do it one cell at a
> time.
> If you get different results, it would be useful to investigate how
> exactly
> they are different.
>
> For continuous elements, the projection is not exact in general, and
> whatever
> you compute locally on one cell has to be reconciled with what you compute
> on
> neighboring cells. The global project is one way to do this reconciliation.
>
> Best
>   W.
>
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/e91a3cb7-d193-b5e4-850f-900ca6a5dea3%40colostate.edu
> .
>

-- 
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/CAM50jEt2KaUEPkA78k4mfCY8xv%3Dzs_SZs8bh19o4b5G9VkXFmQ%40mail.gmail.com.

Re: [deal.II] SparseDirectUMFPACK - member function "initialize" takes very long

[Simon Wiesheier]

I failed to have a look at the wright function. So of couse this answers my
first question.

I measure the times in release mode.

-Simon

Am Mi., 2. Juni 2021 um 16:33 Uhr schrieb Praveen C :

> The initialize function does the LU factorization which is the most
> expensive part
>
>
> https://www.dealii.org/developer/doxygen/deal.II/classSparseDirectUMFPACK.html#af23a74dd9cac006c9d87c18e5d638076
>
> The solution part should be faster than initialize.
>
> Did you measure the time in “release” mode or “debug” ?
>
> best
> praveen
>
> On 02-Jun-2021, at 7:56 PM, Simon  wrote:
>
> Dear all,
>
> I use the SparseDirectUMFPACK solver for doing some post-processing steps.
> I also measure some cpu_times in my program in order to compare different
> methods regarding efficiency. For my question only the following three
> lines of code from my program are relevant:
> ---
> timer.enter_section("Initialize Direct Solver");
> sparse_direct_solver.initialize(global_matrix);
> timer.leave_subsection("Initialize Direct Solver");
> ---
> So I only measure the time which is neeed for executing the member
> function "initialize":
> ->Solving a scalar problem with *dim=2* and *263425 DoFs* the above call
> takes about 1.3 seconds, whereas solving the linear system takes only 0.8
> seconds.
> ->Solving the same scalar problem with *dim=3* and only *72369 DoFs* the
> above call takes about 4.4 seconds and solving the linear system about 3.2
> seconds.
>
> -My first question is what this function actually does since it takes
> longer than assembling and solving the system together?
> The documentation says: "This function does nothing. It is only here to
> provide a interface consistent with other sparse direct solvers."
> By having a look in the function body I´ve seen that this function
> acutally does nothing, i.e. the function body is empty.
>
> -My second question is why for my dim=3 case the function call takes much
> longer as in the dim=2 case? In the latter I have much more DoFs.
> My guess was that the number of DoFs is the essential quantity but this is
> obviously not the case.
>
> Any input would be greatly appreciated!
>
> Best
> Simon
>
>
> --
> 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/5B31FFA2-C446-41CD-A43E-9196A73B8269%40gmail.com
> <https://groups.google.com/d/msgid/dealii/5B31FFA2-C446-41CD-A43E-9196A73B8269%40gmail.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEsKUQitCRov1RuM8%2BEoRMvJjCCsZ8Q%2BGXYDpQe-ooauZw%40mail.gmail.com.

Re: [deal.II] Re: SparseDirectUMFPACK - member function "initialize" takes very long

[Simon Wiesheier]

Thanks for all your answers!
It was my lack of mathematical background regarding solving linear systems
which brought me to that question.

-Simon

Am Mi., 2. Juni 2021 um 17:09 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:

> On 6/2/21 8:59 AM, Abbas wrote:
> >
> > 2) Node numbering plays a role. I am not sure to what extent and I am
> not sure
> > if it is something that is done automatically by direct solvers. Step-2
> with
> > it's Cuthill_McKee
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FnamespaceDoFRenumbering.html%23ab938a690bf4e2adff191fe969b0f21d3&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ca666a43196f746c74d5008d925d6fd84%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637582429206627171%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=FM0KPcJ7TmnahoGrINbqZBr1hL77L8tcSf9ri4DKsHI%3D&reserved=0>
>  is
>
> > something you might want to have a look at .
>
> Yes, all direct solvers I know of renumber automatically. In the case of
> UMFPACK, it uses the minimum-degree renumbering internally -- i.e., there
> is
> no need to do that on the user side.
>
> Best
>   W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/d1a526af-5945-48e4-b3a8-b9ba498b84f7%40colostate.edu
> .
>

-- 
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/CAM50jEtvH0_JOzqA_vu4c9imDDajwAvvBVSy1HYcPgWg1wQFpw%40mail.gmail.com.

Re: [deal.II] Re: SparseDirectUMFPACK - member function "initialize" takes very long

[Simon Wiesheier]

One more question to the renumbering schemes: I called
DoFRenumbering::Cuthill_McKee(dof_handler) before calling the initialize
function of the direct solver.

So this isn´t really necessary when using UMFPACK because internally a
different renumbering scheme is called anyway?

-Simon

Am Mi., 2. Juni 2021 um 17:16 Uhr schrieb Simon Wiesheier <
simon.wieshe...@gmail.com>:

> Thanks for all your answers!
> It was my lack of mathematical background regarding solving linear systems
> which brought me to that question.
>
> -Simon
>
> Am Mi., 2. Juni 2021 um 17:09 Uhr schrieb Wolfgang Bangerth <
> bange...@colostate.edu>:
>
>> On 6/2/21 8:59 AM, Abbas wrote:
>> >
>> > 2) Node numbering plays a role. I am not sure to what extent and I am
>> not sure
>> > if it is something that is done automatically by direct solvers. Step-2
>> with
>> > it's Cuthill_McKee
>> > <
>> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FnamespaceDoFRenumbering.html%23ab938a690bf4e2adff191fe969b0f21d3&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ca666a43196f746c74d5008d925d6fd84%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637582429206627171%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=FM0KPcJ7TmnahoGrINbqZBr1hL77L8tcSf9ri4DKsHI%3D&reserved=0>
>>  is
>>
>> > something you might want to have a look at .
>>
>> Yes, all direct solvers I know of renumber automatically. In the case of
>> UMFPACK, it uses the minimum-degree renumbering internally -- i.e., there
>> is
>> no need to do that on the user side.
>>
>> Best
>>   W.
>>
>> --
>> 
>> Wolfgang Bangerth  email: bange...@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.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/d1a526af-5945-48e4-b3a8-b9ba498b84f7%40colostate.edu
>> .
>>
>

-- 
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/CAM50jEtuqJuvrWSj9gYLLXQ5zgd-%2B6bGCD2hC9u%2BBBMnmuA02g%40mail.gmail.com.

Re: [deal.II] Re:

[Simon Wiesheier]

Hello,

I am dealing with a smilar question as Shahin, i.e. computing global
L2-projections and approximations of mass-matrices with a lumped
integration. So I thought it makes sense to continue in this chat.

I produced two fields via L2-projection:
-In the first case I computed the mass matrix as a matrix, i.e.
(phi_i,phi_j).
-In the second case I approximated the mass matrix  with a lumped
integration. Instead of storing a mass matrix I computed a "mass-vector" by
summing up all off-diagonal elements in each row of the corresponding
mass-matrix and adding it do the corresponding diagonal. From the start I
directly used a vector instead of a matrix.

My distribution of the field (stresses) in terms of critical peaks,... is
of course the same.
However I am a little bit surprised that the "lumped nodal values" result
in lower stresses at most places. I actually expected the other way round,
because:
When thinking on a simple bar-element with 2 nodes, a lumped integration
concentrates the total mass of the bar to the nodes so that the inertia
becomes bigger compared with other Quadrature formulas.

So intuitively speaking, what is the effect of the lumped integration? From
a mathematical viewpoint there are no more couplings between the DoFs. But
does this lead to lower nodal values, higher nodal values or is this
context dependent?

Best
Simon


Am Di., 25. Mai 2021 um 11:52 Uhr schrieb Shahin Heydari <
sh1990.heyd...@gmail.com>:

> Thank you so much Sir, appreciate your help.
>
> Regards
> Shahin
>
> On Mon, May 24, 2021 at 5:46 PM Wolfgang Bangerth 
> wrote:
>
>> On 5/21/21 1:52 PM, Shahin Heydari wrote:
>> > m_i = sigma (m_ij) = int(phi_i dx)
>>
>> Right. You compute the integral as a sum over the integrals on each cell.
>> So
>> when you do
>>
>>for (cell = ...)
>>  for (i=...)
>>{
>>  lump = 0;
>>  for (q=...)
>>lump += ...
>>
>>  lump_cell_matrix(i,i) = lump;
>>}
>>
>> then this is correct if you then add lump_cell_matrix(i,i) to your global
>> matrix.
>>
>> You mention in your original question that you are unsure about what to
>> do
>> with Q2 (or higher) elements. This is a difficult question, but there is
>> a
>> good amount of literature on the topic. Just search for "mass lumping
>> higher
>> order elements".
>>
>> Best
>>   W.
>>
>> --
>> 
>> Wolfgang Bangerth  email: bange...@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.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/1b0908b3-e43c-2d6a-b2a0-09a7a51b4636%40colostate.edu
>> .
>>
> --
> 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/CALUbCFwyACZFLqF2ZEXztJ3H4b2iEhcVGYJ90beAMNPnUVz6aQ%40mail.gmail.com
> <https://groups.google.com/d/msgid/dealii/CALUbCFwyACZFLqF2ZEXztJ3H4b2iEhcVGYJ90beAMNPnUVz6aQ%40mail.gmail.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEtwjTj-BF58cS%3Dn%3DL0E52tFNfbrWYTzspJRnaFRuszfbw%40mail.gmail.com.

Re: [deal.II] Re:

[Simon Wiesheier]

" I'm sure someone has proved this, but my intuition is that for every
elliptic
equation without an advection term, you will get the same convergence rate
as
for the Laplace equation as long as the coefficient in the elliptic
operator
is nice enoug."

I'll will search in the literature to get deeper information, so thanks for
the hint!
And the nonlinear equations of elasticity are actually elliptic equations?

One last question regarding the convergence rate:
As I previously said in this chat, I am projecting the dim*dim coefficients
of the stress tensor, which are only availabe at the qps, to the nodes
using a global L2-projection and for comparison a superconvergent approach
(superconvergent patch recovery introduced by Zienkiewicz and Zhu).
I will then also determine the convergence rate of the two approaches. As
the name suggests, the superconvergent approach should give me a higher
convergence rate
*as expected. *
-Is the *expected rate* the rate which the laplace equation (assuming that
this holds for elasticity as just discussed) yields for the gradients of
the solution, i.e. ||grad (u-u_h)||_L2 <= ...h^p ?

For me it's not yet fully clear how convergence rate of a projection method
(L2, SPR) is linked to ||grad (u-u_h)||_L2 <= ...h^p .
Basically the input for the projection methods are the gradients (stresses)
at the qps. Of course these values will be changed by the method( apply
projection -> get_function_values), but is this "change" totally
independent of the rate p from ||grad (u-u_h)||_L2 <= ...h^p ?
For me this seems to be the case since in the literature people compare the
rate of  a L2-projection with the rate of a SPR-projection (or other
methods).

Thanks again for helping me!

Best
Simon

Am Mo., 7. Juni 2021 um 17:47 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:

> On 6/7/21 8:12 AM, Simon wrote:
> > Yes, actually both schemes converge to the same solution.
> > Maybe I should do some postprocessing with a visualiuation in order to
> figure
> > out what the lumped integration does with the values at the qps.
>
> Then the difference you observe is simply numerical error.
>
>
> > I want to determine the convergence rate for the elasticity equations
> with the
> > PDE: div(stress-tensor) = 0.
> > -Since I determine the convergence rate, the before mentioned *constant*
> in
> > the inequality cancels out, i.e. I should get the same convergence rate
> as for
> > the laplace equation. Is that right?
>
> Correct. The *rate* is independent of the factor.
>
>
> > -I actually solve the nonlinear elasticity equations (hyperelasticity);
> The
> > finite element spaces, test functions,... are the same as for linear
> > elasticity, but the stresses depend now nonlinear on the gradient of the
> > displacements. My question is if the nonlinearity changes the error
> estimates?
> > Let´s assume I determine the convergence rate in the L2-norm for linear
> > elasticity, which is in the best case 2 for for p=1. Should I also get 2
> for
> > the nonlinear elasticity equations?
>
> I'm sure someone has proved this, but my intuition is that for every
> elliptic
> equation without an advection term, you will get the same convergence rate
> as
> for the Laplace equation as long as the coefficient in the elliptic
> operator
> is nice enough. Nice enough here will mean that your stress-strain
> relationship does not degenerate, i.e.,
>||sigma|| / ||eps||
> does not go to zero or infinity, and is a smooth function of some sort.
>
> Best
>   W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/e751f64b-5cc3-d823-0855-5e51d6dd4d46%40colostate.edu
> .
>

-- 
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/CAM50jEuZUUR0U65WcYWXYynGkt8GfnBDOvHEZh96jNR8FXPTNg%40mail.gmail.com.

Re: [deal.II] Calculating normal vectors to deformed surface with pushforward operation

[Simon Wiesheier]

Hello,

my question fits well in this chat, I am also calculating traction vectors
but using linear theory.

I apply a pressure load on some faces of my body. During assembly of course
this is considered as a neumann boundary condition. The direction of my
pressure load is equal to the normal vectors of that  faces.
When my computation is finished I am recomputing the traction vector on the
same faces where I prescribed it during assembly, i.e. on the qps of the
faces I am computing the following:
cauchy_stress_tensor*Normal_vector,
where the cauchy_stress_tensor is computed by the consitutive relation for
linear theory. So I don't have to distinguish between material and spatial
configuartion by making small strain assumptions.

My problem is that the prescribed traction in the assembly routine is not
the same as computed afterwards, there are deviations up to 5-10%.
But they should actually be the same, should they?

It would be great to get that confirmed, so that I have a starting point
for finding possible error sources.

Best
Simon



Am Di., 11. Mai 2021 um 12:10 Uhr schrieb Alex Cumberworth <
alexandercumberwo...@gmail.com>:

>
> Hi Jean-Paul,
>
> Thanks for your response. The spatial normal vectors are now correct. It
> seems my lack of grounding in differential geometry was the culprit here.
>
> However, the integrated forces are still not as I expect, but I think this
> is a separate issue, so I will create a new thread.
>
> Best,
> Alex
> On Tuesday, May 4, 2021 at 6:58:33 p.m. UTC+2 Jean-Paul Pelteret wrote:
>
>> Hi Alex,
>>
>> Well, the one thing that I can clearly identify as being problematic is
>> that the normal vector does not transform with the deformation gradient
>> (directly), but rather with its cofactor. This is what Nanson’s formula for
>> the transformation of surface areas tell us. We actually have a function
>> that does this transformation, namely
>> Physics::Transformations::nansons_formula()
>> <https://dealii.org/current/doxygen/deal.II/namespacePhysics_1_1Transformations.html#aa82925b742c3708f625a48a7abc440bc>
>>  .
>> Maybe you could try to use that and see if you get the result that you’re
>> looking for.
>>
>> Best,
>> Jean-Paul
>>
>> On 4. May 2021, at 12:36, Alex Cumberworth 
>> wrote:
>>
>> Hello,
>>
>> I would like to calculate traction/stress vectors for the deformed
>> configuration. I have solved for a non-linear elasticity problem and am
>> dealing with the Green-Lagrange strain tensor and the second Piola-Kirchoff
>> stress tensor. I thought the simplest way to get the traction vectors would
>> be to first calculate them in the material reference frame (by calculating
>> the stress tensor and applying it to normal vectors in the material
>> reference), and then using a pushforward operation (by applying the
>> deformation gradient) to produce the stress vectors in the spatial
>> reference frame (and also dividing by the determinant of the deformation
>> vector to take care of the change in area).
>>
>> However, the results do not appear as expected. To understand what was
>> going on, I first decided to apply the pushforward operation to the normal
>> vectors. I used the following code to carry this out.
>>
>> template 
>> class SpatialNormalVectorPostprocess: public DataPostprocessorVector
>> {
>>   public:
>> SpatialNormalVectorPostprocess();
>> virtual void evaluate_vector_field(
>> const DataPostprocessorInputs::Vector& input_data,
>> std::vector>& computed_quantities) const;
>> };
>>
>> template 
>> SpatialNormalVectorPostprocess::SpatialNormalVectorPostprocess():
>> DataPostprocessorVector(
>> "spatial_normal",
>> update_gradients | update_normal_vectors) {}
>>
>> template 
>> void SpatialNormalVectorPostprocess::evaluate_vector_field(
>> const DataPostprocessorInputs::Vector& input_data,
>> std::vector>& computed_quantities) const {
>>
>> for (unsigned int i {0}; i != input_data.normals.size(); i++) {
>> Tensor<2, dim, double> grad_u {};
>> Tensor<2, dim, double> identity_tensor {};
>> for (unsigned int d {0}; d != dim; d++) {
>> grad_u[d] = input_data.solution_gradients[i][d];
>> identity_tensor[d][d] = 1;
>> }
>> const Tensor<2, dim, double> deformation_grad {
>> grad_u + identity_tensor};
>> const Tensor<1, dim, double> material_normal_vector {
>> 

[deal.II] Re: Question about deal.II flexibility

[Konrad Simon]

Dear Abdo,

I do think that Deal.II can help you with solving your problem. If you are 
just starting with Deal.II have a look at the tutorials 
. It is always 
a good idea to start with some code basis that other people have shared and 
that was thoroughly tested. As to your future problem, please have a look 
at the code gallery 
 as well. 
There are codes that other users shared dealing with visco-plastic 
materials, in particular (but not only Jean-Paul).

Cheers,
Konrad

On Saturday, June 19, 2021 at 10:23:37 PM UTC+2 aael...@ncsu.edu wrote:

> Hi
> I am Abdo, Ph.D. student at NC state university and was looking for a FE 
> package library that is flexible and can be used in modeling and analyzing 
> 2D and 3D elastic wave equation propagation within an incompressible 
> medium. Also, at some point, I might need to add the viscoelastic behavior 
> to the model. So, I am not sure if deal.II is flexible enough to add, 
> modify some feature in it or not such that I can do what I have mentioned 
> earlier. So, Can you please give me some pointers or advice regarding 
> this?, your help is much appreciated.
>
> Best regards.
> Abdo.
>

-- 
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/481afb9d-fa92-4eae-b369-0881cbd7d70fn%40googlegroups.com.

Re: [deal.II] Re: writing values in the solution vector for constrained dofs

[Simon Wiesheier]

Hi Bruno,

everything worked as expected, so constraints.distribute() does what it
should do :-)

Best
Simon

Am Mi., 23. Juni 2021 um 23:28 Uhr schrieb Bruno Turcksin <
bruno.turck...@gmail.com>:

> Simon,
>
> Yes, you should be able to skip the constrained hanging nodes. Like you
> said the value will be overwritten when you call distribute(). Let us know
> if you get something unexpected but that should work.
>
> Best,
>
> Bruno
>
> On Wednesday, June 23, 2021 at 5:00:06 AM UTC-4 Simon wrote:
>
>> Dear all,
>>
>> I am dealing with adaptivity, so hanging nodes are present in my mesh
>> which are handled by a constraints object like in the following:
>> AffineConstraints...constraints;
>> make_hanging_node_constraints...();
>> constraints.close();
>>
>> My solution vector is filled via a local approach. Basically I loop over
>> all active vertices, use cell->vertex_dof_index() to access the dofs living
>> on them and write appropriate values in the global solution vector. So I
>> don't use the before mentioned constraints object to make a local-global
>> distribution or something similar as it is done in most tuturial programs.
>> After my loop is finished I call constraints.distribute().
>>
>> So by looping over all active vertices I also visit the constrained
>> hanging nodes and write some values for these dofs in the solution vector.
>>
>> My question is if I can just skip the hanging nodes in my loop because
>> the call "constraints.distribute()" will overwrite these values anyway?
>> This is at least my understanding of what constraints.distribute() does,
>> i.e. computing constraind dof values from unconstrained ones. But I am not
>> quite sure if problems can come up by modifyiing the constrained dofs
>> beforehand.
>>
>> Best
>> Simon
>>
> --
> 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/6288fd1c-fa80-4798-ac7b-49a904f1ab40n%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/6288fd1c-fa80-4798-ac7b-49a904f1ab40n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEsbgWAHFPc7QXvqd0ePjgmsXM-xAW1tW9jOCoktp_9peQ%40mail.gmail.com.

Re: [deal.II] Re: std::set> : missing operator(s)

[Simon Wiesheier]

Hi Peter,

Thanks a lot!

It worked for me like this (described in
https://www.cplusplus.com/reference/set/set/set/ ):
bool(*fn_pt)(Point,Point) = comp_points; //comp_points is my
function as I declared it in the question
std::set< Point, bool(*)(Point,Point) > vertices(fn_pt);

Just for a better understanding, the following approach from your linked
website did not work:
std::set< Point, decltype(&comp_points) > vertices(&comp_points);

There are four errors:
-decltype cannot resolve address of overloaded function
-template argument 2 is invalid
-request for member ‘insert’ in ‘vertices’, which is of non-class type ‘int’
-cannot resolve overloaded function ‘comp_points’ based on conversion to
type ‘int’

As far as I understood, decltype(&comp_points) returns '
bool(*)(Point,Point ' ,i.e it should be exactly the same as my
first version, isn't it?

Best
Simon


Am So., 27. Juni 2021 um 12:25 Uhr schrieb Peter Munch <
peterrmue...@gmail.com>:

> Hi Simon,
>
> you don't need to modify the point class. You can simply pas to the
> constructor of the set a comparator. See the first answer in
> https://stackoverflow.com/questions/2620862/using-custom-stdset-comparator
> .
>
> Hope this helps,
> PM
>
> On Sunday, 27 June 2021 at 12:15:40 UTC+2 Simon wrote:
>
>> Dear all,
>>
>> I am trying to create a set of Points, i.e std::set>.
>> I get a quite extensive error message when trying to insert a point in
>> the set but for my understanding of the message, the only problem is the
>> missing operator< for two points.
>> (Well I could make use of a std::vector instead but a std::set fits much
>> better for my issue since (i) I need to insert points at arbitrary
>> positions quite often and second (ii) wanna make sure that a given Point is
>> only once in the container.)
>>
>> But irrespective of what container I use, I would like to supplement the
>> Point Class by adding the 'operator<' like in the following:
>>
>> template
>> bool Point::operator< (const Point& p)
>> {
>> //...some Asserts
>> return this.norm()> }
>>
>> Of course I can not write this piece of code since there is no
>> declaration for a member function 'operator<' in the Point Class.
>> For that reasons my next idea was to modify the file <
>> deal.II/base/point.h
>> <https://www.dealii.org/current/doxygen/deal.II/base_2point_8h_source.html>>
>> by simply adding the missing operator< declaration and the corresponding
>> definition. So I did this but without effect, probably because all the
>> files of the dealii-library are already compiled and executing my program
>> via 'make run' did not see my changes in the Point Class at all. (I also
>> added a dummy command, std::cout<<"...Test..."<> appear on my screen either.)
>>
>> I am only familiar with basic cmake- and make-commands.
>> My question is if it is possible to recompile one specific file of the
>> library again?
>> So if this were possible I would appreciate getting some guidance for
>> implementing this.
>> I would also be interested in writing a patch for that issue.
>>
>> Best
>> Simon
>>
> --
> 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/34faa5ae-9bed-4a66-ad9f-cb2bb59d7ba1n%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/34faa5ae-9bed-4a66-ad9f-cb2bb59d7ba1n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEsthHLQf8nb_qUL3ziAebTh4V%2B6p2wQ_C-aUPyKHvMk9g%40mail.gmail.com.

Re: [deal.II] Re: solution of linear system using LAPACKFullMatrix - problems with singular matrix

[Simon Wiesheier]

Thank you for your reply!

I scaled the x and y coordinates so that they are both between -1 and 1.
Actually this brought the condition numbers down to 1e-04.
I am just a little bit surprised that my results are nearly the same as
those obtained with reciprocal conditon numbers smaller than my machine
accuracy, i.e. about 1e-18.

Would you say that one can continue work with reciprocal condition numbers
of 1e-04 or is this still too big?

But anyway, glad to see that my program is no longer aborting because LU
fails. :-)

Best
Simon

Bruno Turcksin  schrieb am Mi., 30. Juni 2021,
14:48:

> Simon,
>
> Unfortunately, there is no such thing as a preconditioner for direct
> solver. The best way to fix your problem is to try to scale x and y so that
> they have the same order of magnitude. Another way is to try SVD
> <https://dealii.org/current/doxygen/deal.II/classLAPACKFullMatrix.html#aa855f1fb4ccc0d3d263382d8545bbc6a>.
> You may need to check that you get an accurate solution for the finest mesh.
>
> Best,
>
> Bruno
>
> On Tuesday, June 29, 2021 at 6:23:32 AM UTC-4 Simon wrote:
>
>> Dear all,
>>
>> I have to deal with a large number of very small linear systems A*a=b; I
>> need to solve this system for each vertex of my triangulation.
>> The matrix A is purely determinand by dyadic products of certain
>> quadrature points.
>> Example: Point qp =[x,y]; Vector qp_ = [1, x, y, xy, x^2,
>> y^2]^T
>> -> A = A + qp_*qp_^T
>> The number of qps is restricted to the qps of those cells which are in
>> the patch. For instance there are 16 qps when using Q2 elements in 2d with
>> reduced integration (QGauss(2)).
>>
>> First I decided to actually invert the matrix A,
>> i.e. A.invert(A),
>> but the computation of the inverse introduced too much numerical errors.
>> Hence I decided to switch to LAPACKFullMatrix, computed LU factorization
>> and finally solved the system via the .solve() member functions.
>>
>> Switching to LAPACK made my results much better but I still reach a point
>> (fine meshes) where the LU decomposition fails. For nearly each vertex the
>> reciprocal condition number reaches the machine accuracy.
>> I guess the problem is that fine meshes lead to large differences in the
>> x- and y-coordinates of  the qps (Point dummy = [x,y] = [10.0 ,0.01])
>> so that some entries in the A matrix are close to zero and others in turn
>> become very large. IMHO the LU decomposition becomes more demanding if the
>> entries in A differ by several orders. Here is a typical example for A
>> (lowest order e-01 , highest order e+02):
>>
>> 1.600e+01  6.099e+01  3.798e+00  2.327e+02  1.203e+00  1.449e+01
>>  6.099e+01  2.327e+02  1.449e+01  8.884e+02  4.594e+00  5.532e+01
>>  3.798e+00  1.449e+01  1.203e+00  5.532e+01  4.292e-01  4.594e+00
>>  2.327e+02  8.884e+02  5.532e+01  3.395e+03  1.755e+01  2.114e+02
>>  1.203e+00  4.594e+00  4.292e-01  1.755e+01  1.632e-01  1.640e+00
>>  1.449e+01  5.532e+01  4.594e+00  2.114e+02  1.640e+00  1.755e+01
>>
>> My question is if I can give some "massage" to A before calling the LU
>> decomposition? What I mean is something like a preconditioner which are
>> often used when talking about iterative solvers. But when checking the
>> documentation I couldn't find certain functionalities for the Matrix
>> classes.
>> For me it would be great to know if the A matrix could be "rescued"
>> somehow by certain solution techniques, preconditioning or whatever.
>> Otherwise I have to think about the way of constructing A itself.
>>
>> Thanks for giving me feedback!
>>
>> Best
>> Simon
>>
> --
> 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 a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/IEhvdDouaBc/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/a67c534a-c7da-472a-99fa-72d02990c10en%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/a67c534a-c7da-472a-99fa-72d02990c10en%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEvTSwjD3ZaXx20n7UaAbyPxjyxGPYGQ8fE27-q0D0dE0A%40mail.gmail.com.

Re: [deal.II] triangulation.signals.post_refinement :connection of a member function

[Simon Wiesheier]

Thank you very much! With that hint it works as expected.
Is it also possible to get a reference (pointer) to the Solid object within
the constructor of 'AnotherClass', so without passing a solid reference to
the latters constructor? 'this' would in this case refer to the instance of
'AnotherClass'.

The class 'AnotherClass' needs access to objects of type DoFHandler
and CellDataStorage. As for the DoFHandler the copy
instructor, i.e something like 'dof_handler(dof_handler_from_solid)' will
delete 'dof_handler_from_solid'. It actually doesn't delete it but throws
an exception as it is also described in the documentation.
I checked a few source code files of the library how they manage the
storage of large objects, for instance the GridTools::Cache class
stores a triangulation like this:
SmartPointer>  tria;
So I did this in the same way for the DoFHandler and the CellDataStorage
object:
SmartPointer> dof_handler;
...
Is this a efficient way to do that or is there a more sophisticated way to
store the member variables? I asked myself why dealii did not just use the
SmartPointers of the c++ standard library, namely std::shared_ptr<...>,
std::unique_ptr<...>,  I guess most of them also make sure that they
don't become dangling pointers.

Thanks again!



Am Di., 6. Juli 2021 um 05:20 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:

> On 7/5/21 2:56 PM, Simon wrote:
> >
> > So the constructor of Solid initializes "pointer_M" which points to an
> object
> > of "AnotherClass". The variable "tria_signal" in the line below is the
> one
> > which causes the call of the member function
> > &*AnotherClass::mark_for_update.
> > *
> > The above code compiles and actually calls *mark_for_update*  from
> > *AnotherClass* .
> >
> > But what I'd like to achieve is that the signal triggers a member
> function of
> > 'Solid', i.e. something like &*Solid::call_this_function. *So a
> naive
> > approach is to simply replace &*AnotherClass::mark_for_update*  by
> > &*Solid::call_this_function* . (Both functions take no parameters
> and
> > return void).
> > Unfortunately this change produces two error messages.
> > (i) error: pointer to member type ‘void (Solid<2>::)()’ incompatible
> with
> > object type ‘AnotherClass<2>’
> > (ii) error: return-statement with a value, in function returning 'void'
> > [-fpermissive]
> >
> > My understanding of the error message is that only member functions of
> > 'AnotherClass' can be connected to the signal.
> > Is this true?
> > And can I modify the parameters of the connect function somehow to
> actually
> > connect a member function of Solid?
>
> No, but you need to provide a 'this' pointer for the object you want to
> work
> on. In the first case, this is pointer_M.get(). For a member function of
> the
> current object, it would probably be 'this'.
>
> Best
>   W.
>
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/bd80797b-029d-a63e-1124-52a69237bab5%40colostate.edu
> .
>

-- 
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/CAM50jEuzLS3n2STSQJ_AdiMs9rpEd81O8USRZ4aEA7fhSaxnbw%40mail.gmail.com.

Re: [deal.II] triangulation.signals.post_refinement :connection of a member function

[Simon Wiesheier]

I will look that up :-)

Just one last question:
In the constructed a class I store a DoFHandler and Triangulation object as
member variables like this:
SmartPointer> dof_handler;
SmartPointer> triangulation;

I created the instance at the same time when constructing the Solid
instance (pM = boost::make_shared<...>)
After a few loadsteps my solid triangulation changes and the DoFHandler as
well (standard adaptivity).
The instance of the other class where I stored the just mentioned
SmartPointers is "aware" of that, i.e. I don't update the SmartPointers
after adaptivity.
Of course the (smart)pointer stores the address of the triangulation, but
is it not possible that the the call "triangulation.execute_refinement()"
makes the triangulation object too big so that it has to be reallocated at
another adress (->address changes)?
So is that what I observe something like undefined behaviour or do calls
like "dof_handler.distribute_dofs(fe)" or
"triangulation.execute_refinement()" indeed don't change the address of the
underlying triangulation (dof_handler) in memory?

Thank you very much again!



Am Mi., 7. Juli 2021 um 06:35 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:

> On 7/6/21 2:45 AM, Simon Wiesheier wrote:
> >
> > Is it also possible to get a reference (pointer) to the Solid object
> within
> > the constructor of 'AnotherClass', so without passing a solid reference
> to the
> > latters constructor? 'this' would in this case refer to the instance of
> > 'AnotherClass'.
> >
> > The class 'AnotherClass' needs access to objects of type DoFHandler
> and
> > CellDataStorage. As for the DoFHandler the copy
> instructor,
> > i.e something like 'dof_handler(dof_handler_from_solid)' will delete
> > 'dof_handler_from_solid'. It actually doesn't delete it but throws an
> > exception as it is also described in the documentation.
> > I checked a few source code files of the library how they manage the
> storage
> > of large objects, for instance the GridTools::Cache class stores a
> > triangulation like this:
> > SmartPointer>  tria;
> > So I did this in the same way for the DoFHandler and the CellDataStorage
> object:
> > SmartPointer> dof_handler;
> > ...
> > Is this a efficient way to do that or is there a more sophisticated way
> to
> > store the member variables? I asked myself why dealii did not just use
> the
> > SmartPointers of the c++ standard library, namely std::shared_ptr<...>,
> > std::unique_ptr<...>,  I guess most of them also make sure that they
> don't
> > become dangling pointers.
>
> All of this is possible with minimal pain by using "lambda functions" and
> "captured variables". You might want to read up on them in a recent C++
> book :-)
>
> Cheers
>   W.
>
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/0bad1924-dbba-c826-d91e-ed56c17767a6%40colostate.edu
> .
>

-- 
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/CAM50jEuuDrsx0i6SFOw6fTUw3ByTVXLRvEyeJEiP2PX9iivMWQ%40mail.gmail.com.

[deal.II] Re: profiling and parallel performance of deal.II user codes

[simon...@gmail.com]

Hi,

I've used the scalasca tool suite for profiling:

https://www.scalasca.org/scalasca/about/about.html

which use score-p for measuring and the cube gui for visualizing the 
results:

https://www.vi-hps.org/projects/score-p
https://www.scalasca.org/scalasca/software/cube-4.x/download.html

I think it is really good. However, it is somewhat difficult to get started 
with. If you are interested, this user guide is a good place to start:

https://www.scalasca.org/scalasca/software/scalasca-2.x/documentation.html

Best,
Simon

-- 
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/e8ebed6a-4bd9-418c-9fb9-751f9b2babc7n%40googlegroups.com.

Re: [deal.II] finite element with shape functions having delta_ij property at quadrature points

[Simon Wiesheier]

"You've already found this: Both the FE_Q and FE_DGQArbitraryNodes classes
have
constructors that create shape functions based on this information."

The latter is the element I was looking for. The problem  with the
former is that the first and last quadrature point has to be located at
zero and one, respectively - a condition which is not satisfied in my case.
That said, FE_DGQArbitraryNodes is the right element for me, because I do
not use this element to produce a global (dis)continuous field but only do
a local approximation on each cell individually. In particular, I do not
need to create a DoFHandler for this element or do calls like
constraints.distribute_local_to_global().

Although FE_DGQArbitraryNodes works fine for me, I asked myself why the
constructor of FE_Q wants a quadrature object whos first and last qp is at
zero and one, respectively? One can not use this constructor for the family
of QGauss objects, just to mention one instance.

Best
Simon

Wolfgang Bangerth  schrieb am So., 15. Aug. 2021,
20:00:

>
> > I am solving a problem in 2d using FE_Q(2) elements and a gauss
> quadrature
> > rule with (fe.degree +1) quadrature points in each co-ordinate
> direction, that
> > is, I have in total nine quadrature points. My question pertains to the
> following:
> > At each cell, I need to approximate a field whose sampling (support)
> points
> > are the quadrature points belonging to reduced integration, i.e, there
> are
> > four quadrature points in my case and the four (shape) functions
> approximating
> > my field should be designed as follows:
> > N_j (xi_k) = delta_{jk} ,
> > where xi_k are the coordinates of the quadrature points. So I need four
> > (shape) functions, each of which is one at one of the four quadrature
> points
> > and zero at the three others.
>
> You've already found this: Both the FE_Q and FE_DGQArbitraryNodes classes
> have
> constructors that create shape functions based on this information.
>
> > That said, my ansatz is given by (the coefficents a(xi) are of course
> known)
> > f(xi) = a(xi_1) N_1(xi) + a(xi_2) N_2(xi) + a(xi_3) N_3(xi) + a(xi_4)
> N(xi)
> > and I need to evaluate the function f(xi) at the *nine* quadrature
> points.
> >
> > What is the best way to do set up the FEValues object?
> > I have seen that there is a constructor for the FE_Q element which takes
> a
> > Quadratute<1> object. I guess this would help me to define the (shape)
> > functions pertaining to the field f(xi), but I think I can not evaluate
> this
> > field at the nine quadrature points, because (i) their local coordinates
> are
> > of course different in the new FEValues object and (ii) second, I would
> have
> > to insert negative local coordinates for a set of them.
> > Maybe I do not even need a FEValues object for my purpuse. As I said, I
> only
> > need to do the approximation f(xi) and evaluate it at the nine
> quadrature points.
>
> The construction of a finite element field and its evaluation at
> quadrature
> points are two different things. Let's say you use one of the classes
> above to
> create a finite element with the delta-property you seek. Then you create
> a
> DoFHandler with it that describes a finite element field on the entire
> mesh.
> To evaluate it at certain points, you'd just create a FEValues object as
> always, which allows you to obtain the values of shape functions and of
> the
> finite element field that results from a solution vector, at the
> quadratrure
> points of interest.
>
> Best
>   W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/3432bf6b-3048-10df-8ee3-9751bf74c847%40colostate.edu
> .
>

-- 
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/CAM50jEvFrZa53-abnmvWhY9o29AY1rGyfXizwyr8ishBmBqUrA%40mail.gmail.com.

[deal.II] Re: shift-invert transformation problem with slepc

[simon...@gmail.com]

Hi,

I think you need to set both which eigen pair and the target value. For 
example, if you look for the eigenvalue closest in magnitude to 0, you need 
to set the following

const double target_eigenvalue = 0;
eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
eigensolver.set_target_eigenvalue(target_eigenvalue);

Depending on slepc version it might be that you can only choose 
EPS_TARGET_MAGNITUDE. See question 8 here:

https://slepc.upv.es/documentation/faq.htm

Best,
Simon

On Tuesday, August 24, 2021 at 3:39:29 AM UTC+2 anton.i...@gmail.com wrote:

> Dear All,
>
> I am trying to use a shift-invert transformation in step-36. I did minimal 
> changes in the solve function adding the following code that describes what 
> spectral transformation I need to do:
>
> SolverControl solver_control (dof_handler.n_dofs(), 1e-9);
> SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control);
> eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);
> eigensolver.set_problem_type (EPS_GHEP);
>
> double eigen_shift = 1.0e-4;
> SLEPcWrappers::TransformationShiftInvert::AdditionalData 
> additional_data(eigen_shift);
> SLEPcWrappers::TransformationShiftInvert shift (MPI_COMM_WORLD, 
> additional_data);
> eigensolver.set_transformation (shift);
>
> eigensolver.solve(stiffness_matrix,
>   mass_matrix,
>   eigenvalues,
>   eigenfunctions,
>   eigenfunctions.size());
>
> and it gives me the following error: 
>
> [0]PETSC ERROR: Shift-and-invert requires a target 'which' (see 
> EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 
> -eps_target_magnitude
>
> Any help to resolve this error is greatly appreciated. I also attached the 
> full code of modified step-36.
>
> Thank you,
> Anton.
>

-- 
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/ada0420b-83b5-4af0-a77c-d71c9af63dd6n%40googlegroups.com.

[deal.II] Re: Munmap_chunk(): invalid pointer

[simon...@gmail.com]

Hi,
"Invalid pointer" suggest that this is some bug related to memory. One way 
to find such errors is to turn on the adress sanitizer. If you use gcc or 
clang, this can be done by configuring cmake this:

cmake -DCMAKE_BUILD_TYPE="Debug" -DCMAKE_CXX_FLAGS="-fsanitize=address" 
-DCMAKE_EXE_LINKER_FLAGS="-fsanitize=address" ..

and then recompiling. If you now run your program, it will crash with 
something like this:

=
==29575==ERROR: AddressSanitizer: heap-buffer-overflow on address 
0x60207a34 at pc 0x55a3f14c9cb4 bp 0x7ffd346a5f40 sp 0x7ffd346a5f30
WRITE of size 4 at 0x60207a34 thread T0
#0 0x55a3f14c9cb3 in Poro::Poroelasticity<2>::assemble_system() 
/home/simon/workspace/sandbox/example.cc:249
#1 0x55a3f14c110f in Poro::Poroelasticity<2>::run() 
/home/simon/workspace/sandbox/example.cc:394
#2 0x55a3f14b517a in main /home/simon/workspace/sandbox/example.cc:411
#3 0x7fea8d0a80b2 in __libc_start_main 
(/lib/x86_64-linux-gnu/libc.so.6+0x270b2)
#4 0x55a3f14b4f5d in _start 
(/home/simon/workspace/sandbox/build/example+0x76f5d)


At the top entry #0, you can see that the code crashed at line 249,  which 
corresponds to the line:

freq[j]=j+1;

Here, freq had length 1 and you by mistake tried to access the j=1 entry.

Best,
Simon

On Monday, September 13, 2021 at 10:49:41 AM UTC+2 masia...@gmail.com wrote:

> Hello everyone, 
> I am trying to run the following script, but getting the error: 
> 'munmap_chunk(): invalid pointer'. The eclipse debugger doesn't help me 
> realize, where exactly the error occurs. 
> In case you have an idea, how to correct it, please, let me know. I would 
> be very gratefull!
>
> Kind regards,
> Mariia 
>
>
>
>
> #include 
> #include 
> #include 
>
> #include 
> #include 
> #include 
> #include 
> #include 
>
> #include 
> #include 
>
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
>
> #include 
> #include 
>
> //#include 
>
> #include 
> #include 
> #include 
> #include 
>
> #include 
>
> namespace Poroelasticity
> {
>   using namespace dealii;
>
>   template 
>   SymmetricTensor<4, dim, std::complex> 
> get_stress_strain_tensor(const std::complex lambda,
>const std::complex mu)
>   {
>   const std::complex nothing(0.0,0.0);
>   SymmetricTensor<4, dim, std::complex> stress_strain_tensor;
>   for (unsigned int i=0; i< dim; ++i)
>   for (unsigned int j=0; j< dim; ++j)
>   for (unsigned int k=0; k< dim; ++k)
>   for (unsigned int l=0; l< dim; ++l)
>   stress_strain_tensor[i][j][k][l]= (((i==k) && (j==l) ? mu : nothing) +
>   ((i==l) && (j==k) ? mu : nothing) +
> ((i==j) && (k==l) ? lambda : nothing));
>   return stress_strain_tensor;
>   }
>
>   template 
>   inline SymmetricTensor<2, dim> get_strain(const FEValues &fe_values,
>   const unsigned int shape_func,
> const unsigned int q_point)
>   {
>   SymmetricTensor<2, dim> strain_tensor;
>   for (unsigned int i=0; i< dim; ++i)
>   strain_tensor[i][i] = fe_values.shape_grad_component(shape_func, 
> q_point, i)[i];
>
>   for (unsigned int i=0; i< dim; ++i)
>   for (unsigned int j=i+1; j< dim; ++j)
>   strain_tensor[i][j] =
>   (fe_values.shape_grad_component(shape_func, q_point, i)[j] +
>fe_values.shape_grad_component(shape_func, q_point, j)[i]) /
>2;
>   return strain_tensor;
>   }
>
>
>   template 
>   class Poroelasticity
>   {
>   public:
> Poroelasticity(const unsigned int degree);
> void run();
>
>   private:
> void make_grid_and_dofs();
> void assemble_system();
> void solve();
> void output_results() const;
>
> const unsigned int degree;
>
> Triangulation triangulation;
> FESystem  fe;
> DoFHandlerdof_handler;
>
> BlockSparsityPattern  sparsity_pattern;
> BlockSparseMatrix> system_matrix;
> BlockVector> solution;
> BlockVector> system_rhs;
>
> static const SymmetricTensor<4, dim> stress_strain_tensor;
>
> static constexpr unsigned int start_freq=0;
> static constexpr unsigned int end_freq=2; //[Hz]
>   };
>
> template
>   class RightHandSide : public Function
>   {
>   public:
>   RightHandSide()
>   : Function(1)
>   {}
>
>   virtual double value(const Point &p,
>  const unsigned int component = 0) const override;
>   
>   };
>
>   template 
>   double RightHandSide::value

[deal.II] Re: Munmap_chunk(): invalid pointer

[simon...@gmail.com]

Hi,
Good that you solved you problem.

Address sanitizer is great for finding this type of bugs, but note that you 
don't want to use it when you run the actual simulation. 
In general, memory usage increases (by a  factor ~2 to ~3) and the code 
runs slower (by a factor ~2 to ~5).

Best,
Simon

On Monday, September 13, 2021 at 3:03:04 PM UTC+2 simon...@gmail.com wrote:

> Hi,
> "Invalid pointer" suggest that this is some bug related to memory. One way 
> to find such errors is to turn on the adress sanitizer. If you use gcc or 
> clang, this can be done by configuring cmake this:
>
> cmake -DCMAKE_BUILD_TYPE="Debug" -DCMAKE_CXX_FLAGS="-fsanitize=address" 
> -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=address" ..
>
> and then recompiling. If you now run your program, it will crash with 
> something like this:
>
> =
> ==29575==ERROR: AddressSanitizer: heap-buffer-overflow on address 
> 0x60207a34 at pc 0x55a3f14c9cb4 bp 0x7ffd346a5f40 sp 0x7ffd346a5f30
> WRITE of size 4 at 0x60207a34 thread T0
> #0 0x55a3f14c9cb3 in Poro::Poroelasticity<2>::assemble_system() 
> /home/simon/workspace/sandbox/example.cc:249
> #1 0x55a3f14c110f in Poro::Poroelasticity<2>::run() 
> /home/simon/workspace/sandbox/example.cc:394
> #2 0x55a3f14b517a in main /home/simon/workspace/sandbox/example.cc:411
> #3 0x7fea8d0a80b2 in __libc_start_main 
> (/lib/x86_64-linux-gnu/libc.so.6+0x270b2)
> #4 0x55a3f14b4f5d in _start 
> (/home/simon/workspace/sandbox/build/example+0x76f5d)
>
>
> At the top entry #0, you can see that the code crashed at line 249,  which 
> corresponds to the line:
>
> freq[j]=j+1;
>
> Here, freq had length 1 and you by mistake tried to access the j=1 entry.
>
> Best,
> Simon
>
> On Monday, September 13, 2021 at 10:49:41 AM UTC+2 masia...@gmail.com 
> wrote:
>
>> Hello everyone, 
>> I am trying to run the following script, but getting the error: 
>> 'munmap_chunk(): invalid pointer'. The eclipse debugger doesn't help me 
>> realize, where exactly the error occurs. 
>> In case you have an idea, how to correct it, please, let me know. I would 
>> be very gratefull!
>>
>> Kind regards,
>> Mariia 
>>
>>
>>
>>
>> #include 
>> #include 
>> #include 
>>
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>>
>> #include 
>> #include 
>>
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>>
>> #include 
>> #include 
>>
>> //#include 
>>
>> #include 
>> #include 
>> #include 
>> #include 
>>
>> #include 
>>
>> namespace Poroelasticity
>> {
>>   using namespace dealii;
>>
>>   template 
>>   SymmetricTensor<4, dim, std::complex> 
>> get_stress_strain_tensor(const std::complex lambda,
>>const std::complex mu)
>>   {
>>   const std::complex nothing(0.0,0.0);
>>   SymmetricTensor<4, dim, std::complex> stress_strain_tensor;
>>   for (unsigned int i=0; i< dim; ++i)
>>   for (unsigned int j=0; j< dim; ++j)
>>   for (unsigned int k=0; k< dim; ++k)
>>   for (unsigned int l=0; l< dim; ++l)
>>   stress_strain_tensor[i][j][k][l]= (((i==k) && (j==l) ? mu : nothing) +
>>   ((i==l) && (j==k) ? mu : nothing) +
>> ((i==j) && (k==l) ? lambda : nothing));
>>   return stress_strain_tensor;
>>   }
>>
>>   template 
>>   inline SymmetricTensor<2, dim> get_strain(const FEValues 
>> &fe_values,
>>   const unsigned int shape_func,
>> const unsigned int q_point)
>>   {
>>   SymmetricTensor<2, dim> strain_tensor;
>>   for (unsigned int i=0; i< dim; ++i)
>>   strain_tensor[i][i] = fe_values.shape_grad_component(shape_func, 
>> q_point, i)[i];
>>
>>   for (unsigned int i=0; i< dim; ++i)
>>   for (unsigned int j=i+1; j< dim; ++j)
>>   strain_tensor[i][j] =
>>   (fe_values.shape_grad_component(shape_func, q_point, i)[j] +
>>fe_values.shape_grad_component(shape_func, q_point, j)[i]) /
>>2;
>>   return strain_tensor;
>>   }
>>
>>
>>   template 
>>   class Poroelasticity
>>   {
>>   public:
>> Poroelasticity(const unsigned int degree);
>> void run();
>>
>>   private:
>>

[deal.II] Elasticity Preconditioner

[Konrad Simon]

Dear all, 

A student and me are trying to deal with (parallel, static) linearized 
elasticity similar to step-8.
However, the problem is slightly different since the Lame parameters are 
functions instead of constants. 
Can anyone recommend good solvers & preconditioners for the resulting 
system? Our boundary conditions fix all degrees of freedom in the kernel.

Best,
Konrad

-- 
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/eba0cc0f-6028-4fde-9cee-307e09f8d2ban%40googlegroups.com.

[deal.II] Triangulation from lower-dimensional object

[Konrad Simon]

Dear all,

I was recently looking for a solution to the following problem:

I have a triangulation, e.g., of a unit cube. I need to solve a 
lower-dimensional (FEM-like) problem on all faces with BCs on edges and for 
each face I need to solve a problem on all edges with BCs at vertices.

Can I somehow extract a Triangulation form a Triangulation by 
providing, for example, the boundary id on the Triangulation level?

Best,
Konrad

-- 
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/ae0eae00-a572-4036-b20e-31f45dd70892n%40googlegroups.com.

[deal.II] Re: range-based loop over cells

[simon...@gmail.com]

You will iterate from begin() to end() on the object after the colon. On 
Triangulation this is the cells on level 0.

Best,
Simon

On Thursday, March 24, 2022 at 1:57:29 PM UTC+1 Niklas Fehn wrote:

> thanks. I think it is nevertheless possible to write "cell : 
> triangulation", and I wondered what this means exactly?
>
> bruno.t...@gmail.com schrieb am Donnerstag, 24. März 2022 um 13:35:41 
> UTC+1:
>
>> Niklas,
>>
>> If you want cell_iterators, you write for(auto cell : 
>> triangulation.cell_iterators 
>> ()) If you want active_cell_iterators, you write for(auto cell : 
>> triangulation.active_cell_iterators ()) The documentation is here 
>> <https://dealii.org/current/doxygen/deal.II/group__CPP11.html#gaef378969994082255fbc64366511a7d1>
>> .
>>
>> Best,
>>
>> Bruno
>> On Thursday, March 24, 2022 at 8:21:42 AM UTC-4 Niklas Fehn wrote:
>>
>>> Dear dealii,
>>>
>>> when writing
>>>
>>> for(auto cell : triangulation){} 
>>>
>>> which type of cell iterator is used (cell_iterators(), 
>>> active_cell_iterators())? and where do I find this information in the 
>>> documentation?
>>>
>>> Many thanks and 
>>>
>>> Best
>>> Niklas
>>>
>>

-- 
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/c3238c4b-d16d-4102-a10d-6f66db511ff7n%40googlegroups.com.

[deal.II] Re: range-based loop over cells

[simon...@gmail.com]

No, sorry... All cells. Just starting at level 0. Got confused by this line:

Triangulation 
<https://www.dealii.org/developer/doxygen/deal.II/classTriangulation.html>< 
dim, spacedim >::cell_iterator 
<https://www.dealii.org/developer/doxygen/deal.II/group__Iterators.html#ga997d61ac7cdc2be3ae934b1f7cdb>
 
Triangulation 
<https://www.dealii.org/developer/doxygen/deal.II/classTriangulation.html>< 
dim, spacedim >::begin ( const unsigned int 
<https://www.dealii.org/developer/doxygen/deal.II/classint.html>  *level* = 
0) const

/Simon

On Thursday, March 24, 2022 at 2:10:39 PM UTC+1 simon...@gmail.com wrote:

> You will iterate from begin() to end() on the object after the colon. On 
> Triangulation this is the cells on level 0.
>
> Best,
> Simon
>
> On Thursday, March 24, 2022 at 1:57:29 PM UTC+1 Niklas Fehn wrote:
>
>> thanks. I think it is nevertheless possible to write "cell : 
>> triangulation", and I wondered what this means exactly?
>>
>> bruno.t...@gmail.com schrieb am Donnerstag, 24. März 2022 um 13:35:41 
>> UTC+1:
>>
>>> Niklas,
>>>
>>> If you want cell_iterators, you write for(auto cell : 
>>> triangulation.cell_iterators 
>>> ()) If you want active_cell_iterators, you write for(auto cell : 
>>> triangulation.active_cell_iterators ()) The documentation is here 
>>> <https://dealii.org/current/doxygen/deal.II/group__CPP11.html#gaef378969994082255fbc64366511a7d1>
>>> .
>>>
>>> Best,
>>>
>>> Bruno
>>> On Thursday, March 24, 2022 at 8:21:42 AM UTC-4 Niklas Fehn wrote:
>>>
>>>> Dear dealii,
>>>>
>>>> when writing
>>>>
>>>> for(auto cell : triangulation){} 
>>>>
>>>> which type of cell iterator is used (cell_iterators(), 
>>>> active_cell_iterators())? and where do I find this information in the 
>>>> documentation?
>>>>
>>>> Many thanks and 
>>>>
>>>> Best
>>>> Niklas
>>>>
>>>

-- 
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/8403cace-6ba4-4240-8d83-7d36b1832348n%40googlegroups.com.

[deal.II] MPI: Evaluating several distributed functions on each node

[Konrad Simon]

Dear all,

I am facing a little challenge: Suppose I have several objects myFunction 
of a certain class MyFunction derived from a Function object on 
each process in my MPI universe. The number of MyFunction objects may vary 
from process to process.
I am trying to evaluate each object at certain points p that live on a 
distributed::Triangualtion. The evaluation points are different for 
each MyFunction object. A point may be contained in a cell owned by another 
process whose id is known, so I know where to send it to to ask a 
MyFunction object on that specific process to evaluate it for me and send 
the value back.

Now my question: How do I deal with several MyFunction objects on one 
process (I can identify these by a hash). Is it a good idea to have a class 
with a set of static variables and methods on each node to manage the 
communication and distribute necessary information within the process using 
the specific MyFunction hash? Does that make sense?

Let me stress again, the difficulty here is that I have several objects 
that may ask for values one by one or even in a (shared memory) parallel 
manner that may or may not trigger communication (sort of some_to_some 
style as implemented in Deal.II). 

Thanks in advance :-)

Best,
Konrad

-- 
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/4fb6b294-aa7f-4dc2-b823-e42247bf2302n%40googlegroups.com.

Re: [deal.II] How to use create_mass_matrix tools with FESystem

[Simon Sticko]

Hi,

In the stacktrace it looks like the error does not come from create_mass_matrix 
but from VectorTools::project. It complains that the passed Function has 1 
component when it should have dim components.

If you are projecting your own function, you need to make sure you set the 
number of components in the constructor, e.g. as for the body-force in step-18:

https://github.com/dealii/dealii/blob/39e29ba7e1cb8c7f711df6b89d19a76394b16f6a/examples/step-18/step-18.cc#L552

Best,
Simon



On 10/04/2022 01:51, Matthew Rich wrote:

Hi all,

I am trying to use the helper tools to assemble a system for a vector valued 
problem.

I am looking to have 2 or 3 displacement DoFs.

So I use

FESystem fe;

and initialize with

fe(FE_Q 
<https://www.dealii.org/current/doxygen/deal.II/classFE__Q.html>(1), dim)

So my code compiles, but when I try and run I get the following


An error occurred in line <182> of file 

 in function
     void dealii::VectorTools::internal::project_matrix_free(const dealii::Mapping&, const dealii::DoFHandler&, 
const dealii::AffineConstraints&, const dealii::Quadrature&, const dealii::Function::value_type>&, dealii::LinearAlgebra::distributed::Vector&, bool, const 
dealii::Quadrature<(dim - 1)>&, bool) [with int components = 2; int fe_degree = 1; int dim = 2; Number = double; int spacedim = 2; typename 
dealii::LinearAlgebra::distributed::Vector::value_type = double]
The violated condition was:
     dof.get_fe(0).n_components() == function.n_components
Additional information:
     Dimension 2 not equal to 1.

Stacktrace:
---
#0  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free<2, 1, 2, double, 2>(dealii::Mapping<2, 
2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::LinearAlgebra::distributed::Vector::value_type> const&, 
dealii::LinearAlgebra::distributed::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#1  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free_degree<2, 2, double, 2>(dealii::Mapping<2, 
2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::LinearAlgebra::distributed::Vector::value_type> const&, 
dealii::LinearAlgebra::distributed::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#2  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free_component<2, double, 2>(dealii::Mapping<2, 
2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::LinearAlgebra::distributed::Vector::value_type> const&, 
dealii::LinearAlgebra::distributed::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#3  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free_copy_vector<2, dealii::Vector, 
2>(dealii::Mapping<2, 2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints::value_type> 
const&, dealii::Quadrature<2> const&, dealii::Function<2, dealii::Vector::value_type> const&, dealii::Vector&, bool, 
dealii::Quadrature<(2)-(1)> const&, bool)
#4  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project, 2>(dealii::Mapping<2, 2> 
const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints::value_type> const&, dealii::Quadrature<2> 
const&, dealii::Function<2, dealii::Vector::value_type> const&, dealii::Vector&, bool, dealii::Quadrature<(2)-(1)> 
const&, bool)
#5  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::project<2, dealii::Vector, 2>(dealii::Mapping<2, 2> const&, 
dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints::value_type> const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::Vector::value_type> const&, dealii::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#6  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::project<2, dealii::Vector, 2>(dealii::DoFHandler<2, 2> 
const&, dealii::AffineConstraints::value_type> const&, dealii::Quadrature<2> const&, dealii::Function<2, 
dealii::Vector::value_type> const&, dealii::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#7  ./final_proj: final_proj::WaveEquation<2>::run(

Re: [deal.II] Encountering error in installation of deal.ii

[Simon Sticko]

Hi,

From the error message, it looks like you don't have C++ compiler installed.
If you are using ubuntu, try

sudo apt-get install g++

and rerun cmake.

Best,
Simon

On 08/06/2022 16:32, Abhishek Nath Thakur wrote:

Respected deal.ii users,

I'm trying to install deal.ii 9.2.0 version but keep getting the error(*error 
image attached below*). I'm very new to this so i can really use someone's help 
to find me the way out of this.

Thank you very much in advance


--
The deal.II project is located at http://www.dealii.org/ 
<http://www.dealii.org/>
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en 
<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 
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e5dd7aeb-ade5-4c4a-9774-6ce10535b3b3n%40googlegroups.com
 
<https://groups.google.com/d/msgid/dealii/e5dd7aeb-ade5-4c4a-9774-6ce10535b3b3n%40googlegroups.com?utm_medium=email&utm_source=footer>.


--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6e38f5a1-4032-47c5-e31c-0036353adc35%40gmail.com.

Re: [deal.II] Facing problem while installation deal.ii

[Simon Sticko]

Hi,

Just to elaborate slightly on what Chen said: If your installation stops at 52%, it is NOT because 
of "dependencies not satisfied". This is not an error. Can you post the error you get 
from calling make in the terminal? E.g. copy the output you get from cmake and "make" to 
a txt-file and attach that.

Best,
Simon

On 10/06/2022 10:05, 陈敏 wrote:

Dear Abhishek Nath Thakur

These are many examples showing how to use deal.ii, these are corresponding documents on 
the official website <https://dealii.org/developer/doxygen/deal.II/Tutorial.html>. 
But some of these need external packages 
<https://dealii.org/developer/readme.html#:~:text=bit%2Dindices.-,Optional%20interfaces%20to%20other%20software%20packages,-When%20configuring%20interfacing>
 not bundled in the deal.ii. you need to install them yourself. Finally, you don't need to 
install all of the external packages, only the necessary packages are needed according to 
your applications.

Best
Chen

Abhishek Nath Thakur mailto:abhisheknaththa...@gmail.com>> 于2022年6月10日周五 15:57写道：

Dear deal.ii users,

I tried to install deal.ii 9.2.0 several times but every time i'm getting 
"dependencies not satisfied" (*image attached below*) due to which my 
installation stops at 52% everytime. I'm not sure if its an error or something else. So 
any help is appreciated

Thank you in advance


-- 
The deal.II project is located at http://www.dealii.org/ <http://www.dealii.org/>

For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en 
<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 
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5483bfa8-a269-4503-8624-557e8aa015edn%40googlegroups.com
 
<https://groups.google.com/d/msgid/dealii/5483bfa8-a269-4503-8624-557e8aa015edn%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
The deal.II project is located at http://www.dealii.org/ 
<http://www.dealii.org/>
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en 
<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 
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADr3OdL74ZJdfrwXhe1M%2BB0sZ4%2BwtNBBhV4bGUCqv1vspUvEqg%40mail.gmail.com
 
<https://groups.google.com/d/msgid/dealii/CADr3OdL74ZJdfrwXhe1M%2BB0sZ4%2BwtNBBhV4bGUCqv1vspUvEqg%40mail.gmail.com?utm_medium=email&utm_source=footer>.


--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1bb4cf10-c237-938f-666a-7b87e175d13d%40gmail.com.

Re: [deal.II] Integral over part of domain.

[Simon Sticko]

Hi,

This type of integral can be computed with the tools presented in the Step-85 
tutorial:

https://www.dealii.org/developer/doxygen/deal.II/step_85.html

which you can use if you install the master branch of deal.II (or wait for the 
soon-released 9.4 version).

In particular, take a look at the assemble_system() function in Step-85 and the 
NonMatching::FEValues class:

https://www.dealii.org/developer/doxygen/deal.II/classNonMatching_1_1FEValues.html

Best,
Simon


On 13/06/2022 15:17, Oleg Kmechak wrote:

Hello there,

Facing such a problem.
Want to evaluate integral over circle (see picture below). Integral has such a 
form:
I = Integral[ solution_value(x, y) * my_func(x, y) * dx dy,  over circle].

Currently I am using the simplest way to evaluate such integral - rectangular 
approximation. And most time cost part - is using 
Functions::FEFieldFunction to get solution_value(x, y). Solving FEM is 
faster than evaluation of such integrals.

Also, my_func(x, y) has hidden parameter 'order'. Basically, I am using 
integral to expand field around tip (in the center of circle) into series. And 
also with higher order, accuracy of series parameter (result of integration) is 
poorer.

So maybe once again. How to evaluate fast and accurately integral over such 
circle (not whole domain)?

Also, maybe I can introduce a function which is equal to 1 inside of circle and 
0 - outside. And then integrate over whole domain.
Also, maybe I can adopt The deal.II Library: Integrators (dealii.org) 
<https://www.dealii.org/current/doxygen/deal.II/group__Integrators.html> to fit 
my function but not sure how I can do it.

Best regards,
Oleg Kmechak

Inkedinit_river_with_boundary_cond_LI.jpg

--
The deal.II project is located at http://www.dealii.org/ 
<http://www.dealii.org/>
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en 
<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 
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e02ebf7a-c4f2-4202-8cd2-057900260914n%40googlegroups.com
 
<https://groups.google.com/d/msgid/dealii/e02ebf7a-c4f2-4202-8cd2-057900260914n%40googlegroups.com?utm_medium=email&utm_source=footer>.


--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/79693521-76dc-ece1-e0eb-32fb0bf7c36c%40gmail.com.

Re: [deal.II] How to add degrees of freedom to a specified node

[Simon Sticko]

Hi,

Since several XFEM papers formulate this as a Heaviside enrichment around the 
crack, you might find this class useful:

https://www.dealii.org/current/doxygen/deal.II/classFE__Enriched.html

Relating to Wolfgangs answer, note in particular the section "Using enriched and 
non-enriched FEs together".

Best,
Simon


On 18/06/2022 07:43, Wolfgang Bangerth wrote:


Wang Yuan,


  I am using extended finite element method to do crack expansion. I have a 
question, can I increase the degree of freedom for the nodes designated around 
the crack in Deal?


Not easily. deal.II enumerates degrees of freedom based on the finite element 
used on a cell, so if you want to have extra degrees of freedom, you need to 
describe that by using a different element on that cell (using the 
hp-framework, in the same way as step-46 for example).


Meanwhile, another question is can I assign different numbers of Gaussian 
integral points to different elements? Thank you for your help!


Yes! You can use a hp::QCollection together with hp::FEValues. Or you can do it 
by hand: You create multiple quadrature rules and corresponding FEValues 
objects, and then you choose which one is appropriate when you are on a 
concrete cell.

Best
  W.



--
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/56f2ee5e-9bba-fd46-7499-7b925aa3353d%40gmail.com.

Re: [deal.II] mesh_generator

[Simon Sticko]

Hi,

Did you call

triangulation.execute_coarsening_and_refinement();

after setting the refinement flag?

Best,
Simon

On 23/06/2022 17:22, LY XXXiao wrote:

Hi there,

Has anyone experienced this issue in mesh refinement:

I have already a basic mesh with 120 cells, and I only want to refine one 
element among them, so I set a refine_flag for this element. But the refinement 
does not happen ?? Is there anything else that I missed? I learned step 1 and 
step 6 for meshing problem.

Best regards,
Longying


--
The deal.II project is located at http://www.dealii.org/ 
<http://www.dealii.org/>
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en 
<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 
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0cfddd51-4380-4cd3-81af-9259046b2c8en%40googlegroups.com
 
<https://groups.google.com/d/msgid/dealii/0cfddd51-4380-4cd3-81af-9259046b2c8en%40googlegroups.com?utm_medium=email&utm_source=footer>.


--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/49efb643-4ca8-f2a6-c9a7-c7c20fa12ee8%40gmail.com.

Re: [deal.II] shape function gradient at arbitrary points in the element

[Simon Wiesheier]

Dear Martin,

thanks for pointing out to the FEPointEvaluation class.
However, my intent is to compute the gradients of shape functions in real
coordinates.

" so in that case the recommendation would be to use FEValues, despite that
being terribly inefficient."

What is the approach to achieve this?
My only idea is to create a new FEValues object for each quadrature point
and pass an appropriate Quadrature object. But seems to be very inefficient.
Is there an alternative way in combination with FEValues?

Best,
Simon


Am Di., 5. Juli 2022 um 08:16 Uhr schrieb Martin Kronbichler <
kronbichler.mar...@gmail.com>:

> Dear Simon,
>
> We have the class FEPointEvalation,
> https://dealii.org/developer/doxygen/deal.II/classFEPointEvaluation.html
> , which implements the operation of evaluating a solution at arbitrary
> points (as handed in as an array to points in unit coordinates). To do this
> for positions in real coordinates, you typically need to invert the
> mapping, Mapping::transform_real_to_unit_cell or
> Mapping::transform_points_real_to_unit_cell, see
> https://dealii.org/developer/doxygen/deal.II/classMapping.html , always
> assuming that you have located the correct points. Most tasks can be
> reduced to this setup.
>
> If your intent is not just the solution interpolation but the value of
> shape functions, e.g. because you want to assemble a local matrix,
> FEPointEvaluation is not directly providing that info, so in that case the
> recommendation would be to use FEValues, despite that being terribly
> inefficient. We could possibly integrate such functionality into
> FEPointEvaluation if you have interest in that direction, so let us know
> what your needs are.
>
> Best,
> Martin
> On 05.07.22 15:56, Simon wrote:
>
> Dear all:
>
> I have to compute the gradient with respect to real space co-ordinates of
> the shape functions belonging to a FE_Q element at arbitrary points in the
> element, that is, not only at quadrature points.
>
> What is the way to do this in dealii?
> FEValues provides shape function values, gradients,... only at quadrature
> points.
>
> Thank you,
> Simon
> --
> 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/16afc8fa-bcfe-4a81-9ade-040c18f85596n%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/16afc8fa-bcfe-4a81-9ade-040c18f85596n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/34843cd0-22a6-d4b5-7cb8-f530715a6108%40gmail.com
> <https://groups.google.com/d/msgid/dealii/34843cd0-22a6-d4b5-7cb8-f530715a6108%40gmail.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEteHa8vhryPshNi4SGP3jhm%3DsEj_A8K6wq1TSX-aamKmQ%40mail.gmail.com.

Re: [deal.II] discretization of a analytical function on a second triangulation

[Simon Wiesheier]

Dear Daniel,


yes, I defined a class derived from deal::Function and it seems to work for
me.

Basically, my motivation to call VectorTools::interpolate is to compute the
nodal values (scalar field --> DOFs) in the triangulation based on the
formula implemented in the just derived class.
During the assembly routine, I access some interpolated values at
quadrature points,...

With this in mind, is this an appropriate way?


Best
Simon

Am Di., 5. Juli 2022 um 09:16 Uhr schrieb Daniel Arndt <
d.arndt.m...@gmail.com>:

> Simon,
>
> defining a class derived from dealii::Function sounds sensible.
> Whether calling VectorTools::interpolate makes sense, depends on what you
> want to do with the interpolated function. In may cases, it's enough to
> just evaluate it locally during assembly using functionality like
> FEValues::get_function_values.
>
> Best,
> Daniel
>
> On Mon, Jul 4, 2022 at 10:26 AM Simon  wrote:
>
>> Dear all:
>>
>>
>> I am solving a nonlinear PDE on a Triangulation T_1.
>>
>> I know the analytical representation of a scalar function of two
>> variables (no space co-ordinates, but two invariants of a quantity)
>> and my goal is to find a discretized version of the analytical function
>> (on a second, two-dimensional triangulation T_2).
>> The nodal values of T_2 should be the exact values corresponding to the
>> analytical solution.
>> As for the values between the nodes, I want to start with a bilinear
>> interpolation.
>>
>> The coupling to my original PDE works as follows:
>> For each quadrature point defined on T_1, I compute the two input
>> variables of the above function based on my solution vector, find the
>> corresponding cell of T_2 and want to retrieve the interpolated nodal
>> values.
>>
>> My idea is to use the function
>> VectorTools::interpolate ,
>> but I do not know a suitable Function object to hand over.
>>
>> Is this an appropriate approach at all?
>>
>>
>> Thank you,
>> Simon
>>
>> --
>> 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/14b1d8f6-f902-4414-a5c6-348ae3f5003en%40googlegroups.com
>> <https://groups.google.com/d/msgid/dealii/14b1d8f6-f902-4414-a5c6-348ae3f5003en%40googlegroups.com?utm_medium=email&utm_source=footer>
>> .
>>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/CAOYDWbKNTV%3D8Bipyq_gzbMoLvxdvv-%2BmHLipn27DbOcK-eWL_A%40mail.gmail.com
> <https://groups.google.com/d/msgid/dealii/CAOYDWbKNTV%3D8Bipyq_gzbMoLvxdvv-%2BmHLipn27DbOcK-eWL_A%40mail.gmail.com?utm_medium=email&utm_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAM50jEsXcMsvZzq0-EFgQwWn%3Dr%3D6yBw8yVrMnQ5met_LK0zHew%40mail.gmail.com.

Re: [deal.II] shape function gradient at arbitrary points in the element

[Simon Wiesheier]

Just to make sure that I chose a reasonable approach:

//create FEValues object whenever I need it, that is,
'number_of_cells*number_of_qps_per_cell' times
FEValues<2> fe_values(dof_handler.get_fe(),

Quadrature<2>(my_point_in_ref_coords),
   update_gradients);
fe_values.reinit(cell);

//evaluate, say, the gradient of the second shape function in real coords
at "my_point_in_ref_coords"
Tensor<1,2> shape_gradient_real_second = fe_values.shape_grad(1, 0);

Without measuring the time but based on your intuition - the approach using
FEPointEvaluation would not be much faster, would it?

Best
Simon



Am Do., 7. Juli 2022 um 03:45 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:

> On 7/5/22 09:28, Simon Wiesheier wrote:
> >
> > What is the approach to achieve this?
> > My only idea is to create a new FEValues object for each quadrature
> point and
> > pass an appropriate Quadrature object. But seems to be very inefficient.
>
> This is basically what the VectorTools::point_value() and
> VectorTools::point_gradient() functions do. You might want to read through
> their implementations to see how this is done.
>
> The only other approach is with classes like FEPointEvaluation. If you
> read
> through their implementation, you might be able to understand how they can
> also be used to evaluate just a single shape function. In the end, the
> values
> and gradients of a single shape function correspond to the values and
> gradients of a finite element field that corresponds to a solution vector
> with
> one 1.0 and the rest of the elements set to zero.
>
> Best
>   W.
>
>
> --
> 
> Wolfgang Bangerth  email: bange...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/6527dccb-aebc-f424-5217-1c19a9b8b3c2%40colostate.edu
> .
>

-- 
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/CAM50jEsJh_oDF%3DzZFvnq69rVbcn5GQmBvji8u2oLyeU8SbfK6A%40mail.gmail.com.

Re: [deal.II] automatic differentiation:

[Simon Wiesheier]

Hi Jean-Paul,

thanks for pointing that out.

However, I believe there has been a misunderstanding:
I am working here with a FE discretized energy on a seperate triangulation,
say, 'triangulation_energy', with the coordinates being the first and third
invariant of 'C'; The corresponding nodal energy values are stored in the
vector 'solution_energy' and are computed once based on a analytical energy
function.
So for each quadrature point in my geometry (first triangulation) I am
computing 'C=F^T*F', thence the first and third invariant of 'C' and store
them in a Point<2> object denoted as 'invariant_point'.
Then, I employ
'GridTools::find_active_cell_around_point(triangulation_energy,
invariant_point)' to finally create a second FEValues object with the only
purpose to compute the gradient of the energy ('solution_energy') at
'invariant_point'.
'get_function_gradients(solution_energy, grad_energy_at_ref_point ) gives
me the gradient of 'solution_energy' with respect to the first and third
invariant of 'C' because the real coordinates on the associated
triangulation are the first and third invariant of 'C'.
Now, I have all the quantities to compute the stress tensor as
S_ad = 2*(grad_energy_at_ref_point[0][0]*invariant1_wrt_C_ad
+
grad_energy_at_ref_point[0][1]*invariant3_wrt_C_ad);

My goal is to compute the 4th order tangent \frac{\partial S_ad}{\partial
C_ad}.

All that said, is it not correct to start like this?:
SymmetricTensor<2,dim> C = symmetrize(F^t *F);
ad_helper.register_independent_variable(C, C_dofs);
const SymmetricTensor<2, dim, ADNumberType> C_ad =
ad_helper.get_sensitive_variables(C_dofs);

Best,
Simon



Am Mo., 18. Juli 2022 um 22:46 Uhr schrieb Jean-Paul Pelteret <
jppelte...@gmail.com>:

> Hi Simon,
>
> Unfortunately I don’t have the time at this moment to give you a full
> explanation as to why the logic of your code is wrong, but in essence you
> have the sequence of operations inverted. You need to compute your energy
> based on the “sensitive” DoF values that would come from the reinit’d
> FEValues operation. You cannot compute something like C_ad ahead of time,
> like it appears that you have.
>
> Take a look at these lines in the (unfinished) tutorial that will
> illustrate the steps required to linearise an energy functional. This
> should hopefully help steer you in the right direction.
>
> https://github.com/dealii/dealii/pull/10394/files#diff-fc8b83d1370bfd7eb558ea76175bfc0e8d6305023d54b17ec9cccb0fba9b1e02R1758-R1831
>
> Best,
> Jean-Paul
>
> On 18. Jul 2022, at 19:46, Simon  wrote:
>
> Dear all,
>
> I want to retrieve a tangent at quadrature point level using AD - attached
> is a snippet of my code:
>
> The computation of the dependent variable 'S_ad' in line 36 involves the
> call
> get_function_gradients (solution_energy, grad_energy_at_ref_point),
> see line 35.
> 'solution_energy' is a scalar function which depends on the independent
> variable 'C_ad' as can be seen in the derivations from line 8 on.
> My issue is that I see no way to pass the NumberTypes type
> (=Differentiation::AD::NumberTypes::sacado_dfad_dfad) to the corrosponding
> FEValues object, and, as a consequence, the AD framework is not aware that
> 'solution_energy' is a function of 'C_ad'.
> In particular, 'grad_energy_at_ref_point' will be treated as a constant
> with respect to 'C_ad' in line 36f.
>
> How can I incorporate this dependency in the AD framework?
>
> Thank you,
> Simon
>
> --
> 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/990f7f7f-37a2-4f9d-b6c2-bbf7dda21e8dn%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/990f7f7f-37a2-4f9d-b6c2-bbf7dda21e8dn%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
> 
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to