Re: [deal.II] LinearOperator -- inverse_operator : Payload for LinearAlgebra::distributed::Vector

[Matthias Maier]

Hi Michał,

It would be helpful if you could sent us a small, minimal example that
shows the problem (in particular, what type is "LevelVectorType" and
"SCPreconditoinMGs", "solver").
The example should be as short as possible and doesn't need to be
functional / or do anything at all - we essentially only need to be able
to figure out the types you try to wrap into the linear_operator object.

Best,
Matthias


On Tue, Feb 13, 2018, at 12:47 CST, Michał Wichrowski  
wrote:

> Dear all, 
> I'm trying to use iterative inverse with distributed deal.II vector. From 
> following lines 
>
> const auto S= linear_operator (mg_sc_matrices[level]);
> const auto Shat = 
> inverse_operator(S,solver,SCPreconditionMGs[level]);
>
>  I get compilation errors:
>
>
> /home/mwichro/lib/deal.II/include/deal.II/lac/linear_operator.h:678:1: 
> note: candidate: template Preconditioner, class Range, class Domain> dealii::LinearOperator Range, Payload> dealii::inverse_operator(const 
> dealii::LinearOperator&, Solver&, const 
> Preconditioner&)
>  inverse_operator(const LinearOperator ,
>  ^
> /home/mwichro/lib/deal.II/include/deal.II/lac/linear_operator.h:678:1: 
> note:   template argument deduction/substitution failed:
> /home/mwichro/deal-ii-projects/StokeMatrixFree/StokesMatrixFree.cc:1051:55: 
> note:   mismatched types 
> ‘dealii::LinearAlgebra::distributed::Vector’ and 
> ‘dealii::internal::LinearOperator::EmptyPayload’
> const auto Shat = 
> inverse_operator(S,solver,SCPreconditionMGs[level]);
>^
> /home/mwichro/deal-ii-projects/StokeMatrixFree/StokesMatrixFree.cc:1051:55: 
> note:   ‘const 
> dealii::LinearOperator dealii::LinearAlgebra::distributed::Vector, 
> dealii::internal::LinearOperator::EmptyPayload>’ is not derived from ‘const 
> dealii::LinearOperator dealii::LinearAlgebra::distributed::Vector >’
> CMakeFiles/StokesMatrixFree.dir/build.make:62: recipe for target 
> 'CMakeFiles/StokesMatrixFree.dir/StokesMatrixFree.cc.o' failed
> make[2]: *** [CMakeFiles/StokesMatrixFree.dir/StokesMatrixFree.cc.o] Error 1
>
>
> I think that proper Payload is needed, but I didn't found right one.
>
> Michał 

-- 
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] Wrong sign appears in VectorTools::point_gradient function when dealing with hyperball mesh

[Wolfgang Bangerth]

On 02/16/2018 07:34 PM, Moraad Biagooi wrote:


As you suggested, I tried a 2D mesh. The plots shows that the sign difference
happens at the edge of two neighbor symmetric cells answers. Two gradient
functions may calculate on different cells, so it looks like the signs are 
different.


So it looks to me that there was no problem to begin with.


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

Re: [deal.II] Re: Crack propagation

[feapman]

Dear Prof. Wick,

I have used isotrope formulation for miehe shear loading (without local 
refinement). I ca not see two crack branching which is described in the 
article (A review on phase-field models of brittle fractureand a new fast 
hybrid formulation)

I attach the test results.

Thanks for your answer!

Kind regards,
Yaakov

On Saturday, February 3, 2018 at 7:20:27 AM UTC+1, Thomas Wick wrote:
>
> Hello,
>
> there are two flags in the parameter file that you need to change:
>
> set Decompose stress in rhs   = 0.0
> set Decompose stress in matrix   = 0.0
>
>
> Thomas W.
>
>
> -- 
> ++++
> Prof. Dr. Thomas Wick
> Institut für Angewandte Mathematik (IfAM)
> Leibniz Universität Hannover
> Welfengarten 1
> 30167 Hannover, Germany
>
> Tel.:   +49 511 762 3360
> Email:  thoma...@ifam.uni-hannover.de 
> www:http://www.ifam.uni-hannover.de/wick
> www:http://www.cmap.polytechnique.fr/~wick/
> ++++
> -- 
>
>
>
>
> On 02/03/2018 03:24 AM, fea...@gmail.com  wrote:
>
> Dear Prof. Heister,
>
> I would like to just test  isotrope formulation of phase -field model (no 
> compression/tension modification), how could I modify the codes (in a 
> simple way)?
>
> I am sorry that I am just a beginner of DealII.
>
> Best
>
>
> On Tuesday, December 5, 2017 at 7:01:17 PM UTC+1, Timo Heister wrote: 
>>
>> > If I use your user element, I have to use OPEN MPI? now I have some 
>> issues 
>> > with Open MPI in Deal.ii 
>>
>> What do you mean by "user element"? The example code in question 
>> requires deal.II to be configured with MPI. What vendor you use 
>> (OpenMPI, MPICH, ...) is up to you. 
>>
>> -- 
>> Timo Heister 
>> http://www.math.clemson.edu/~heister/ 
>>
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> --- 
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dealii+un...@googlegroups.com .
> For more options, visit https://groups.google.com/d/optout.
>
>
>

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

[deal.II] GridTools::collect_periodic_faces has problem with meshes generated by Gmsh or Abaqus when running in parallel.

[Hamed Babaei]

Hi All,

I'm having problem when applying periodicity on meshes generated with Gmsh 
or Abaqus. The problem arises when I run in parallel (in serial it works) 
getting the following error:

 void 
dealii::GridTools::match_periodic_face_pairs(std::set >&, std::set >&, int, 
std::vector&, const 
dealii::Tensor<1, typename FaceIterator::AccessorType:: space_dimension>&, 
const dealii::FullMatrix&) [with CellIterator = 
dealii::TriaIterator >; typename 
dealii::identity::type = 
dealii::TriaIterator >; typename 
FaceIterator::AccessorType = dealii::CellAccessor<3, 3>]
The violated condition was:
n_matches == pairs1.size() && pairs2.size() == 0
The name and call sequence of the exception was:
ExcMessage ("Unmatched faces on periodic boundaries")
Additional Information:
Unmatched faces on periodic boundaries 

I suspect that match_periodic_face_pairs function does not have access to 
the whole faces in parallel. 

I am wondering if you have encountered this issue before, or any clue what 
is the reason for this issue.

Thanks in advance.

Hamed

-- 
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] Non smooth solution when the BC coincide with the processors boundaries?

[Giorgos Kourakos]

Thank you for you reply,

The issue with the boundaries was indeed solved when I assigned the 
boundary ids on the ghost cells as well.

As for the other problem with the refinement around the processor borders, 
I was using a relatively complex way to move the triangulation vertices 
each iteration using a solution transfer to interpolate the vertices to the 
refined mesh. 

After a lot of digging I found the method 
communicate_locally_moved_vertices with a note that I should convert the 
triangulation back to its original shape, then apply refinement, and then 
move the vertices to the desired position.

This simple trick works great so far. 
Kudos to whoever came up with that !!

Giorgos

On Thursday, February 15, 2018 at 3:54:03 PM UTC-8, Wolfgang Bangerth wrote:
>
> On 02/15/2018 04:48 PM, Giorgos Kourakos wrote: 
> > 
> > I haven't forgotten the ConstraintMatrix::distribute. The code works 
> > with 1 processor 
> > I do set manually boundary ids but I reset them after every refinement. 
> > Isn't this enough to get the correct boundaries? 
> > I haven't cross out the other two cases yet, although it doesnt seem as 
> > an artifact. I would say that it looks like that the one processor 
> > doesnt know that there is a boundary in this ghost cell. 
> > At the moment I set boundary ids only for the locally owned cells. 
>
> Not good enough. Both the IndexSet you pass to ConstrainedMatrix and 
> VectorTools::interpolate_boundary_values() needs to have correct 
> information about ghost cells. So the index set needs to contain all 
> locally relevant DoFs, and you need to set boundary_ids for all ghost 
> cells as well. 
>
>
> > I'd like to know if this is somehow expected when more than one 
> > processor is used. 
>
> No -- the solution you compute needs to be exactly the same (up to 
> round-off and solver tolerance, of course) whether you compute with one 
> or multiple processors. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

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

Re: [deal.II] Wrong sign appears in VectorTools::point_gradient function when dealing with hyperball mesh

[Wolfgang Bangerth]

On 02/15/2018 07:29 AM, Moraad Biagooi wrote:


I have some problems with the signs of VectorTools::point_gradient. 
Here's an example:


I modified Step-4 of tutorial so that the mesh is a 
GridGenerator::hyperball with half of its boundary_id equal to one.


After solving the problem, a function is called that prints 
VectorTools::point_gradient function versus a simple 1st order gradient 
at some points inside the mesh.


The result is almost the same but there's sign difference at some 
points.


Do you know which function is correct? For example, have you taken a 
look at the solution  at that point and determined visually what the 
gradient is, and then compared it with what the two ways of computing 
things return? (This may be easier to do in 2d.)


There is of course always the possibility of a bug, but it would 
surprise me if there was one in the function that computes the gradient 
at a point.


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

[deal.II] Re: Apply the inhomogeneous Dirichlet boundary conditions on the velocity in the mixed FEM

[jie liu]

Dear Wolfgang,

Thank you for your reply.

This problem is solved by adding 
constraints.distribute (solution);
in 
void MixedLaplaceProblem::solve ()
in my code.

Kind regards,
Jie 

On Thursday, February 15, 2018 at 5:24:44 PM UTC+1, jie liu wrote:
>
>
> Dear all,
>
> I want to solve an 1D Poisson equation with the inhomogeneous boundary 
> condition on the gradient (velocity) using the mixed FEM. The equation reads
>
>
>
> I used step-20 to implement the mixed FEM.
> Since this step cannot deal with the inhomogeneous boundary condtion on 
> the gradient, I added some code of step-22.
> The changes are:
> (1) adding
> system_matrix.clear ();
> dof_handler.distribute_dofs (fe);
> DoFRenumbering::Cuthill_McKee (dof_handler);
> std::vector block_component (dim+1,0);
> block_component[dim] = 1;
> DoFRenumbering::component_wise (dof_handler, block_component);
>
> {
> constraints.clear ();
> FEValuesExtractors::Vector velocities(0);   
> DoFTools::make_hanging_node_constraints (dof_handler,
>  constraints);
> 
> VectorTools::interpolate_boundary_values (dof_handler,
>   1,
>   GradientBoundary(),
>   constraints,
>   fe.component_mask(velocities
> )
>  );
>  
> }
> constraints.close ();
> of  step-22 in
>
>> void MixedLaplaceProblem::make_grid_and_dofs () 
>>
> of step-20;
> (2) adding
> for (unsigned int i=0; i   for (unsigned int j=i+1; j local_matrix(i,j) = local_matrix(j,i);
> cell->get_dof_indices (local_dof_indices);
> 
> constraints.condense (system_matrix, system_rhs);
> of  step-22 in
>   void MixedLaplaceProblem::assemble_system ()
> of step-20.
>
> Pn elements are used for the solution and Pn+1 elements are used for the 
> gradient.
> fe (FE_Q(degree+1), 1,
> FE_Q(degree), 1),
>
> Setting 'a' equal or not equal to zero, it works when 'b' is equal to zero 
> (homogeneous Dirichlet boundary condition on the gradient). 
> Specifically, the convergence order of the solution is one order higher 
> than the element degree of the solution, and the convergence order of the 
> gradient is equal to or one order lower than the element degree of the 
> gradient.
> However, when 'b' is not equal to zero ( inhomogeneous Dirichlet boundary 
> condition on the gradient), it does not work for the gradient, i.e. the 
> convergence order of the solution is still one order higher than the 
> element degree, but the convergence order of the gradient keeps 0.5 for 
> different elements. 
>
> Could anyone explain this? Thank you.
>
> Kind regards,
> Jie
>

-- 
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] Force and PK1

[daschneider95 via deal.II User Group]

> Dear David,
>
> I studied the step- 44 Tutorial Program and read already a lot about how 
> to apply a Neumann contribution. But I have still two questions:
> In my case, I like to apply a force (not a force per unit reference area) 
> as a Neumann contribution on the RHS. But for the Integral, Ni * traction* 
> JxWI need the Piola- Kirchhoff traction.
> My idea would be the following: F= PN dA = T dA   <=> F/dA= T
> I would simply divide the force by the area of the cellsurface. 
> So my first question is: Can I do it like this or is there a simpler way 
> of doing it?
>
>
> It depends. The approach you’ve described here (that is, to specify the 
> total force applied to a surface) is the same as is done in this 
> code-gallery example (see line 1690):
>
> https://github.com/dealii/code-gallery/tree/master/Quasi_static_Finite_strain_Compressible_Elasticity
>
> But this only makes sense if you know the area over which the traction is 
> integrated and if the traction force is uniformly distributed. For a 
> geometrically complex shapes, this may no longer be the best solution, and 
> you should probably consider precomputing the total surface area before 
> defining how the traction is distributed. This point relates to your second 
> question...
>
> And the second question is: How can I access this sort of Information? The 
> area might for example be computed using the norm of a cross product. But I 
> am not sure how to get the two vectors for this operation.
>
>
> You can just use numerical integration for this (which can also take 
> account of the mapping used on curved boundaries). Look at the 
> compute_vol_current() function in step-44, which can be easily adopted to 
> compute surface areas. 
>
> I hope this helps.
>
> Best,
> Jean-Paul
>
>
> First of all thank you very much for the quick replies Lucas and Jean-Paul.

I should have been a little more specific. In my case, I obtain the 
integrated force, which refers to a cellsurface. So, I dont't get the total 
force for my whole surface and the force is only the same for every 
quadrature point of one cellface (uniformly distributed on one cellface). I 
just wanted to make sure that this approach makes sense. For simple 
geometries (no curved boundaries) the function measure() should work fine. 
It might be important to catch also the case of curved edges. I think for 
this case I can integrate the function 1 over the current domain, like 
already mentioned in the descripion of the function measure().

Best regards,
David


 

-- 
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] Force and PK1

[Jean-Paul Pelteret]

Dear David,

> I studied the step- 44 Tutorial Program and read already a lot about how to 
> apply a Neumann contribution. But I have still two questions:
> In my case, I like to apply a force (not a force per unit reference area) as 
> a Neumann contribution on the RHS. But for the Integral, Ni * traction* JxW   
>  I need the Piola- Kirchhoff traction.
> My idea would be the following: F= PN dA = T dA   <=> F/dA= T
> I would simply divide the force by the area of the cellsurface. 
> So my first question is: Can I do it like this or is there a simpler way of 
> doing it?

It depends. The approach you’ve described here (that is, to specify the total 
force applied to a surface) is the same as is done in this code-gallery example 
(see line 1690):
https://github.com/dealii/code-gallery/tree/master/Quasi_static_Finite_strain_Compressible_Elasticity
 


But this only makes sense if you know the area over which the traction is 
integrated and if the traction force is uniformly distributed. For a 
geometrically complex shapes, this may no longer be the best solution, and you 
should probably consider precomputing the total surface area before defining 
how the traction is distributed. This point relates to your second question...

> And the second question is: How can I access this sort of Information? The 
> area might for example be computed using the norm of a cross product. But I 
> am not sure how to get the two vectors for this operation.

You can just use numerical integration for this (which can also take account of 
the mapping used on curved boundaries). Look at the compute_vol_current() 
function in step-44, which can be easily adopted to compute surface areas. 

I hope this helps.

Best,
Jean-Paul


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

[deal.II] Re: Force and PK1

[Lucas Campos]

Dear David,

 
> So my first question is: Can I do it like this or is there a simpler way 
of doing it?
I can't really comment on this. Hopefully someone with more experience than 
me will chip in.
 
> And the second question is: How can I access this sort of Information? 
The area might for example be computed using the norm of a cross product. 
But I am not sure how to get the two vectors for this operation.
You might want to check this 
function: 
http://dealii.org/developer/doxygen/deal.II/classTriaAccessor.html#a07891adb63b1629ef34bf285d1b70af2

I would recommend you read carefully the description of the function before 
using it, mainly for curved domains.

Bests,
Lucas

On Friday, 16 February 2018 08:43:32 UTC+1, daschn...@googlemail.com wrote:
>
> Dear all,
>
> I studied the step- 44 Tutorial Program and read already a lot about how 
> to apply a Neumann contribution. But I have still two questions:
> In my case, I like to apply a force (not a force per unit reference area) 
> as a Neumann contribution on the RHS. But for the Integral, Ni * traction* 
> JxWI need the Piola- Kirchhoff traction.
> My idea would be the following: F= PN dA = T dA   <=> F/dA= T
> I would simply divide the force by the area of the cellsurface. 
> So my first question is: Can I do it like this or is there a simpler way 
> of doing it?
> And the second question is: How can I access this sort of Information? The 
> area might for example be computed using the norm of a cross product. But I 
> am not sure how to get the two vectors for this operation.
>
> Any help is appreciated.
> Thanks a lot with best regards
> David
>

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