Re: [deal.II] Application of the non-homogeneous boundary conditions in Finite Elasticity (Compressible case)

2020-10-11 Thread Jean-Paul Pelteret
Hi Animesh,

Here are some answers to your questions:

> Essentially, I will have to apply the homogeneous boundary condition at the 
> zeroth and first Newton's step, and Non-homogeneous BC only at the zeroth 
> step, as mentioned in the documentation 
> .

> 1. Why would I need to give the homogeneous boundary conditions twice?

No, I’m afraid you’ve misunderstood what’s written there. Here’s the important 
part of the documentation:
“”"
However, since we are dealing with an iterative Newton method, it should be 
noted that any displacement constraints should only be specified at the zeroth 
iteration and subsequently no additional contributions are to be made since the 
constraints are already exactly satisfied.
“””

You only apply (inhomogeneous) displacement constraints at the zeroth Newton 
iteration. Because this tutorial uses an incremental formulation, and the 
constraints are applied to the increment itself, you want to make sure that you 
apply the correct implement only once per timestep. This is because you 
effectively displace your boundary from its position at the previous tilmestep 
to the one at the current tilmestep using this single update. After this update 
is done, you want to make sure that the boundary doesn’t move away from that 
location, so you apply homogeneous constraints on the same displacement 
components for each subsequent Newton step. That is the other part of the logic 
that’s laid out in step-44’s Solid::make_constraints().

> 2. Is it correct the way I am applying the constraints? I am copying the part 
> of the code in red? The major concern that I have is that if I do 
> constraints.clear(), then it will also clear my non-homogeneous constraint 
> applied at the zeroth Newton step. How should I handle this?

Not quite. You do want the inhomogeneities cleared, but you still only want to 
fix the displacement components that you originally applied inhomogeneous 
constraints to. You’ll want to duplicate the logic applied in the other 
branches of the make_constraints() function. So something like this:
> // Apply Non-homogeneous B.C.
> {
>  const int boundary_id = 5;
>   if (apply_dirichlet_bc == true)
> VectorTools::interpolate_boundary_values(dof_handler_ref,
>  boundary_id,
>  ConstantFunction(dim == 
> 3 ? std::vector { -0.0001, 0.0, 0.0 } : std::vector { 
> -0.0001, 0.0 }), constraints, fe.component_mask(x_displacement));
else
> VectorTools::interpolate_boundary_values(dof_handler_ref,
>  boundary_id,
>  ZeroFunction(), 
> constraints, fe.component_mask(x_displacement));
> }


I must admit that this code is overly complicated, and it can be simplified so 
that extending it to the case of inhomogeneous constraints is easier. I’ve 
opened up an issue for this task here 
 and I’ll update the tutorial in 
the coming weeks.

I hope that this helps you in understanding the logic and resolving your issue!

Best,
Jean-Paul


> On 10 Oct 2020, at 19:05, Animesh Rastogi IIT Gandhinagar 
>  wrote:
> 
> Hi All, 
> 
> I am working on hyper-elastic beam buckling problem. I wish to apply 
> non-homogeneous Dirichlet Boundary condition on the right hand side of the 
> beam (to compress the beam). I was reading the step-44 
> 
>  where they have mentioned about how to apply non-homogeneous boundary 
> condition.  Essentially, I will have to apply the homogeneous boundary 
> condition at the zeroth and first Newton's step, and Non-homogeneous BC only 
> at the zeroth step, as mentioned in the documentation 
> .
> 
> I have following two questions regarding this - 
> 
> 1. Why would I need to give the homogeneous boundary conditions twice?
> 2. Is it correct the way I am applying the constraints? I am copying the part 
> of the code in red? The major concern that I have is that if I do 
> constraints.clear(), then it will also clear my non-homogeneous constraint 
> applied at the zeroth Newton step. How should I handle this?
> 
> I have edited the following code from code gallery - 
> 
> https://dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html
>  
> 
> 
> Thanks!
> 
> Animesh
> 
> 
> 
> 
> template 
>   void Solid::make_constraints(const int _nr)
>   {
> std::cout << " CST " << std::flush;
> if (it_nr > 1)
>   return;
> constraints.clear();
> const bool apply_dirichlet_bc 

Re: [deal.II] "locked" vtu outputs

2020-10-11 Thread David Montiel Taboada
Thank you, Timo

I will look into it!

David

On Sat, Oct 10, 2020 at 11:30 AM Timo Heister  wrote:

> Hi David,
>
> > For example, for the initial conditions, the code is supposed to output
> a file with name:
> >
> > solution-000.vtu
> >
> > but then it also creates s file with name
> >
> > solution-000.vtu-1364024289-50151.lock
> >
> > This is what the support from the cluster told us about the issue:
>
> The code you posted uses MPI I/O for parallel output. This lock file
> is created by the operating system and to me it looks like the
> filesystem you are writing to might be an NFS storage that is not set
> up to work well with parallel IO. Generally, NSF is not a great option
> for performance. It might also be that the file system is not
> configured correctly.
>
> I have a few suggestions:
> 1. Use a real parallel file system for IO (most clusters have support
> for a fast parallel file system and a slow home directory on NSF)
> 2. Change the code to not use parallel IO (replace
> write_vtu_in_parallel by a call to  write_vtu_with_pvtu_record() and
> set n_groups = 0). This is slower, especially for large computations
> with 100+ cores. You might need to change a couple of lines of code.
> 3. Talk to the admins to change/update the MPI implementation
> (assuming this actually is a parallel file system)
>
> Best,
> Timo
>
>
> --
> 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+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/CAMRj59E1v%2BN3%3DXfmqLPV3NX_hhNuG_zwoMsejs_2b_6CL1Grng%40mail.gmail.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/CAJdNbL867m4E0GVnq3q7qSVnv%3Dc4BRsD%2BsJN-H5tMug%3D77zh2Q%40mail.gmail.com.


Re: [deal.II] FullMatrix to SparseMatrix

2020-10-11 Thread Nikki Holtzer


I have copied some code below. 


I was hoping to take the outer product which is formed in a FullMatrix, 
called Kronecker_matrix below, multiply it by a coefficient, and add a 
matrix to it called nonlinear_part_full. This entire operation I have tried 
to put into a Sparse matrix, called Nonlinear_part, because I need to add 
this term to linear_part which is Sparse. I have done this simply because 
it was the only way I could make the addition between linear_part (Sparse) 
and nonlinear_part_full (FullMatrix) work all while trying to make the code 
as efficient as possible. 



template 

void TimeStep::prepare()

{ 

  std::cout << "TimeStep::prepare()" << std::endl;

  const auto _data = *p_offline_data;

  linear_part.reinit(offline_data.sparsity_pattern);

  auto _c = offline_data.mass_matrix;

  auto _c = offline_data.stiffness_matrix;

  

  linear_part.copy_from(M_c);

  linear_part.add((1. - theta) * kappa, S_c);

  linear_part_inverse.initialize(linear_part);

  nonlinear_part.reinit(offline_data.sparsity_pattern);

  auto  = offline_data.inner_product_matrix;

  nonlinear_part.copy_from(X);

}


template 

void TimeStep::step(Vector _solution,double new_t)


{

  const auto _data = *p_offline_data;

  const auto  = *offline_data.p_coefficients;

  const auto _constraints = offline_data.affine_constraints;

  auto  = *p_manufactured;


  const auto lam = coefficients.lam;

  std::cout<<"lambda in Time step  = "<> vector_memory;

  typename VectorMemory>::Pointer 
p_new_solution(vector_memory);

  auto _solution = *p_new_solution;

  new_solution = old_solution;

  apply_boundary_values(offline_data, coefficients, new_t, new_solution, 
old_solution);

  

  for (unsigned int m = 0; m < nonlinear_solver_limit; ++m) 
{ Vector residual = M_u * (new_solution - old_solution) 

  + kappa * (1. - theta) * S_u * new_solution +  theta * kappa * 
S_u * old_solution - kappa * (1.- theta) * lam * new_solution.norm_sqr()   
 *  new_solution - kappa * theta * lam * 
old_solution.norm_sqr() * old_solution;


affine_constraints.set_zero(residual);

if (residual.linfty_norm() < nonlinear_solver_tol)

  break;

double old_solution_normsq = -old_solution.norm_sqr();

 const unsigned int k = old_solution.size();

FullMatrixKronecker_matrix(k,k);

for(unsigned int i=0; inonlinear_part_full;

nonlinear_part_full.copy_from(nonlinear_part);

double nonlinear_part_coef = - kappa * (1.- theta) * lam * 
old_solution_normsq;

double Kronecker_matrix_coef = -2.* lam * kappa * (1- theta);

nonlinear_part_full *= nonlinear_part_coef; // -k(1- theta) lambda 
|Psi|^2*X

   Kronecker_matrix *= Kronecker_matrix_coef; //(-2.* lam * kappa * (1- 
theta)* Kronecker_matrix)

   nonlinear_part_full.add(1.,Kronecker_matrix); //-k(1- theta) lambda 
|Psi|^2*X - 2k(1-theta)lambda Kronecker_matrix

   

   SparseMatrix Nonlinear_part;

Nonlinear_part.reinit(offline_data.sparsity_pattern);

Nonlinear_part.copy_from(nonlinear_part_full);

linear_part.add(1., Nonlinear_part);

auto system_matrix = linear_operator(linear_part);



SolverControl solver_control(linear_solver_limit, linear_solver_tol);

SolverGMRES<> solver(solver_control);


auto system_matrix_inverse = inverse_operator(system_matrix, solver, 
linear_part_inverse);

Vector update = system_matrix_inverse * (-1. * residual);


affine_constraints.set_zero(update);



new_solution += update;


}


On Sunday, October 11, 2020 at 5:38:12 PM UTC-4 Wolfgang Bangerth wrote:

>
> Nikki,
> we can continue to speculate, but your description of what is happening is 
> not 
> specific enough to really know. You'll have to show us the code you use to 
> be 
> sure.
>
> I don't understand the connection between the code you show and the outer 
> product. The outer product is a dense matrix, but the sparsity pattern you 
> create is sparse. Were you hoping to put the outer product of two vectors 
> into 
> a sparse matrix?
>
> Best
> W.
>
>
> On 10/10/20 3:27 PM, Nikki Holtzer wrote:
> > My sparsity pattern was created in the following way:
> > 
> > DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), 
> dof_handler.n_dofs());
> > 
> >   DoFTools::make_sparsity_pattern(dof_handler, c_sparsity, 
> > affine_constraints, true);
> > 
> > 
> > which I have used throughout my code with no problems.
> > 
> > 
> > The FullMatrix, Mat1, that I form comes from the outer product of vec1 
> and 
> > vec2 both of size dof_handler.n_dofs.
> > 
> > I create a Sparse Matrix, Mat2 that i initialize with the sparsity 
> pattern 
> > written above. I am not seeing why Mat1 of size dof_handler.n_dofs X 
> > dof_handler.n_dofs does not fit into a sparsity pattern that's defined 
> as 
> > DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), 
> dof_handler.n_dofs());
> > 
> > I am indeed getting an error that there are not enough entries in the 
> > SparseMatrix, Mat2.
> > 
> > Thank you!
> > 
> > 

Re: [deal.II] FullMatrix to SparseMatrix

2020-10-11 Thread Wolfgang Bangerth



Nikki,
we can continue to speculate, but your description of what is happening is not 
specific enough to really know. You'll have to show us the code you use to be 
sure.


I don't understand the connection between the code you show and the outer 
product. The outer product is a dense matrix, but the sparsity pattern you 
create is sparse. Were you hoping to put the outer product of two vectors into 
a sparse matrix?


Best
 W.


On 10/10/20 3:27 PM, Nikki Holtzer wrote:

My sparsity pattern was created in the following way:

DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs());

   DoFTools::make_sparsity_pattern(dof_handler, c_sparsity, 
affine_constraints, true);



which I have used throughout my code with no problems.


The FullMatrix, Mat1, that I form comes from the outer product of vec1 and 
vec2 both of size dof_handler.n_dofs.


I create a Sparse Matrix, Mat2 that i initialize with the sparsity pattern 
written above. I am not seeing why Mat1 of size dof_handler.n_dofs X 
dof_handler.n_dofs does not fit into a sparsity pattern that's defined as 
DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs());


I am indeed getting an error that there are not enough entries in the 
SparseMatrix, Mat2.


Thank you!


On Wednesday, October 7, 2020 at 4:58:45 PM UTC-4 Wolfgang Bangerth wrote:

On 10/7/20 2:51 PM, Nikki Holtzer wrote:
 >
 > When doing so I receive the following error:
 >
 >
 > An error occurred in line <195> of file
 >  in function
 >
 >     dealii::SparseMatrix&
 > dealii::SparseMatrix::operator=(double) [with number = double]
 >
 > The violated condition was:
 >
 >     cols != nullptr
 >
 > Additional information:
 >
 >     (none)

You haven't associated a sparsity pattern with the sparse matrix yet. You'll
have to set that first, and you need to make sure that it has an entry for
each nonzero in the full matrix you are hoping to copy over!

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/56cca082-6ae2-4b51-af0d-85ed95661ac4n%40googlegroups.com 
.



--

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/1c294be4-0d06-0044-0aec-d377c6b54c05%40colostate.edu.


Re: [deal.II] Tools for parallel debugging

2020-10-11 Thread Wolfgang Bangerth

On 10/11/20 3:26 PM, Daniel Arndt wrote:


have a look at 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs 
. 
For me, the "mpiexec -np n gdb ./executable" approach normally works fairly well.




There is also an interactive demonstration of this in lecture 41.25:
https://www.math.colostate.edu/~bangerth/videos.676.41.25.html

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/a1eb1592-a7b1-c471-3693-ce103f5e177b%40colostate.edu.


Re: [deal.II] Tools for parallel debugging

2020-10-11 Thread Daniel Arndt
Konrad,

have a look at
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs.
For me, the "mpiexec -np n gdb ./executable" approach normally works fairly
well.

Best,
Daniel

Am So., 11. Okt. 2020 um 09:56 Uhr schrieb Konrad Simon <
ksimon1...@gmail.com>:

> Hi deal.ii community,
>
> I have a short question about tools (which probably also concerns many
> other people).
>
> I am using eclipse for C++ development and I frequently use the debugger.
> So far everything fine. Now I need to use parallel debugging tools for MPI
> but my eclipse crashes many times and/or is super slow.
>
> What tools are you using? Any recommendations or hints?
>
> 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/d89c4b39-a069-4825-b96b-741665e55095n%40googlegroups.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/CAOYDWbJT%2B6QaCVyOBn3nk408AX0JX%2B_g8A28aZJefZBhKkz75Q%40mail.gmail.com.


[deal.II] Tools for parallel debugging

2020-10-11 Thread Konrad Simon
Hi deal.ii community,

I have a short question about tools (which probably also concerns many 
other people). 

I am using eclipse for C++ development and I frequently use the debugger. 
So far everything fine. Now I need to use parallel debugging tools for MPI 
but my eclipse crashes many times and/or is super slow.

What tools are you using? Any recommendations or hints?

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/d89c4b39-a069-4825-b96b-741665e55095n%40googlegroups.com.


[deal.II] Re: Fluid-Structure interaction

2020-10-11 Thread blais...@gmail.com
Dear Ramiro,
You could also take a look at step-70 by Luca and I. It deplays a very 
prototype-ish fluid-structure interaction framework using immersed 
boundaries. I am sure Luca has progressed a lot on this issue since then. 
You can also look at software derived from Deal.II. Our group develops a 
full Navier-Stokes solver (https://github.com/lethe-cfd/lethe) based on 
deal.II. Extending it to a fluid-structure software could be done 
relatively simply if you want to use a non-monolithic coupling.
Best
Bruno


On Wednesday, October 7, 2020 at 10:29:02 a.m. UTC-4 ramrebol wrote:

> Hi everyone.
>
> I'm looking for a solver to solve a fluid-structure problem, like a heart 
> valve: the valve is an elastic material, and the blood could be considered 
> as an incompressible fluid (evolutionary incompressible Navier--Stokes 
> equations).
>
> In step-46 I see a Fluid-Structure interaction implemented in dealii
>
> https://www.dealii.org/current/doxygen/deal.II/step_46.html
>
> but I'm not sure if it is appropiate to a more realistic problem.
>
> I know that the first thing is to know what is the numerical method that 
> we want implement, but I'm not sure what is the appropiated numerical 
> method to solve this strong coupled problem.
>
> My goal is to implement a model of a valve oppening and close due blood 
> pressure (a cilindrical tube with an elastic obstacle in the middle). Do 
> you think dealii could be an appropiate code to model this? Or is a too big 
> problem?
>
> The other possibility is to use commercial software (COMSOL or ANSIS), but 
> due we don't know what they do to solve the problem, I could prefer to 
> program "my own" solver, for example with dealii.
>
> I suppose there must be many papers proposing algorithms for these types 
> of coupled problems, but I have not found anyone explaining the 
> mathematical aspects behind it. So I am very confused about how to 
> implement this coupled problem.
>
> I don't know if is too much to ask here,  if so, I apologize in advance.
>
> In concrete, my questions are the following:
>
> - what is the numerical algorithm to solve this coupled problem.
> - and, is it possible to solve this with dealii? Or is a too big problem.
>
>
>
>
> -- 
> Ramiro
>

-- 
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/c5a55124-992b-4bec-8135-92bdd60e7c85n%40googlegroups.com.