Re: [deal.II] Step-3 Tutorials : Trying to use bessel function as a boundary dirichlet value

[Pushkar Pandit]

Well I did resolved the above issue by changing the center to a point that
lies on the boundary.

On Fri, Jan 1, 2021 at 10:52 AM pushkar...@gmail.com <
pushkarpandi...@gmail.com> wrote:

> Dear deal.II Community,
>
> In my attempt to learn this library  I was working with step-3 tutorials
> where I wish to change my dirichlet boundary function to a bessel function
> instead of zero function but I ended up with a issue in visualization where
> I am not able to see any variation  on the mesh but I am able to run the
> program successfully.I have added the following code :
>
>  Point<2> center(1.5,1.5);
>  Functions::Bessel1<2> bessel(1,2.5,center);
> VectorTools::interpolate_boundary_values (dof_handler,
> 0,
> bessel,
> boundary_values);
>
> Although I may be wrong conceptually at various places and I look forward
> to any help in this context as my goal with this exercise  to just change
> the homogeneous boundary condition to inhomogeneous one.
>
> Thankfully
> Pushkar
>
>
> --
> 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/8e0029e0-1c64-4d7d-a524-c948682743e5n%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/CAKWKLQJqPjyu0F5HBRQkoT4bnx%2Bfh-p0k1w4sexTkrWZtDOG3Q%40mail.gmail.com.

[deal.II] Step-3 Tutorials : Trying to use bessel function as a boundary dirichlet value

[pushkar...@gmail.com]

Dear deal.II Community,

In my attempt to learn this library  I was working with step-3 tutorials 
where I wish to change my dirichlet boundary function to a bessel function 
instead of zero function but I ended up with a issue in visualization where 
I am not able to see any variation  on the mesh but I am able to run the 
program successfully.I have added the following code : 

 Point<2> center(1.5,1.5);
 Functions::Bessel1<2> bessel(1,2.5,center);
VectorTools::interpolate_boundary_values (dof_handler,
0,
bessel,
boundary_values);

Although I may be wrong conceptually at various places and I look forward 
to any help in this context as my goal with this exercise  to just change 
the homogeneous boundary condition to inhomogeneous one.

Thankfully
Pushkar


-- 
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/8e0029e0-1c64-4d7d-a524-c948682743e5n%40googlegroups.com.

Re: [deal.II] Re: Parallel distributed hp solution transfer with FE_nothing

[Marc Fehling]

Kaushik,

in addition to what I just wrote, your example from above has revealed a 
bug in the `p::d::SolutionTransfer` class that Wolfgang and I were 
discussing in the course of this chatlog. Thank you very much for this! We 
are working on a solution for this issue.

I would encourage you to use the `p::d::CellDataTransfer` class for your 
use case as described in the last message.

Marc

On Thursday, December 31, 2020 at 6:02:00 PM UTC-7 Marc Fehling wrote:

> Hi Kaushik,
>
> Yes, this is possible by changing a cell from FE_Nothing to FE_Q using 
> p-refinement.
>
> You can do this with the method described in #11132 
> : Imitate what 
> p::d::SolutionTransfer is doing with the more versatile tool 
> p::d::CellDataTransfer and consider the following:
>
>- Prepare a data container like `Vector>` where the 
>outer layer represents each cell in the current mesh, and the inner layer 
>corresponds to the dof values inside each cell.
>- Prepare data for the updated grid on the old grid. 
>   - On already activated cells, store dof values with 
>   `cell->get_interpolated_dof_values()`.
>   - On all other cells, store an empty container.
>   - Register your data container for and execute coarsening and 
>refinement.
>- Interpolate the old solution on the new mesh.
>- Initialize your new solution vector with invalid values 
>   `std::numeric_limits::infinity`.
>   - On previously activated cells, extract the stored data with 
>   `cell->set_dof_values_by_interpolation()`.
>   - Skip all other cells which only have an empty container.
>- On non-ghosted solution vectors, call 
>`compress(VectorOperation::min)` to get correct values on ghost cells.
>
> This leaves you with a correctly interpolated solution on the new mesh, 
> where all newly activated dofs have the value `infinity`.
>
> You can now (and not earlier!!!) assign the values you have in mind for 
> the newly activated dofs. You may want to exchange data on ghost cells once 
> more with `GridTools::exchange_cell_data_to_ghosts()`.
>
> This is a fairly new feature in the library for a very specific use case, 
> so there is no dedicated class for transferring solutions across finite 
> element activation yet. If you successfully manage to make this work, would 
> you be willing to turn this into a class for the deal.II library?
>
> Marc
> On Wednesday, December 30, 2020 at 8:22:59 AM UTC-7 k.d...@gmail.com 
> wrote:
>
>> Hi all,
>> Thank you for your reply. 
>> Let me explain what I am trying to do and why. 
>> I want to solve a transient heat transfer problem of the additive 
>> manufacturing (AM) process. In AM processes, metal powder is deposited in 
>> layers, and then a laser source scans each layer and melts and bonds the 
>> powder to the layer underneath it. To simulate this layer by layer process, 
>> I want to start with a grid that covers all the layers, but initially, only 
>> the bottom-most layer is active and all other layers are inactive, and then 
>> as time progresses, I want to activate one layer at a time. I read on this 
>> mailing list that cell "birth" or "activation" can be done by changing a 
>> cell from FE_Nothing to FE_Q using p-refinement. I was trying to keep all 
>> cells of the grid initially to FE_Nothing except the bottom-most layer. And 
>> then convert one layer at a time to FE_Q. My questions are:
>> 1. Does this make sense? 
>> 2. I have to do two other things for this to work. (a) When a cell 
>> becomes FE_Q from FE_Nothing, dofs that are activating for the 1st time, I 
>> need to apply a non-zero initial value to those dofs. This is to simulation 
>> the metal powder deposited at a specified temperature,. e.g. room 
>> temperature. (b) the dofs, that were shared between a FE_Q and FE_Nothing 
>> cells before the p-refinement and now shared between FE_Q and FE_Nothing 
>> cells after refinement, should retrain their values from before the 
>> refinement. 
>>
>> Another way to simulation this process would be to use a cell "awaking" 
>> process, instead of the cell "birth". I keep call cells FE_Q but apply a 
>> very low diffusivity to the cells of the layers that are not yet "awake". 
>> But this way, I have to solve a larger system in all time steps. I was 
>> hoping to save some computation time, by only forming a system consists of 
>> cells that are in the "active" layers only. 
>>
>> Please let me if this makes sense? Is there any other method in deal.ii 
>> that can simulation such a process? 
>> Thank you very much and happy holidays.
>> Kaushik 
>>
>>
>> On Tue, Dec 29, 2020 at 12:26 PM Wolfgang Bangerth  
>> wrote:
>>
>>> On 12/28/20 5:11 PM, Marc Fehling wrote:
>>> > 
>>> > In case a FE_Nothing has been configured to dominate, the solution 
>>> should be 
>>> > continuous on the interface if I understood correctly, i.e., zero on 
>>> the face. 
>>> > I will write a few tests to see if this is a

Re: [deal.II] Re: Parallel distributed hp solution transfer with FE_nothing

[Marc Fehling]

Hi Kaushik,

Yes, this is possible by changing a cell from FE_Nothing to FE_Q using 
p-refinement.

You can do this with the method described in #11132 
: Imitate what 
p::d::SolutionTransfer is doing with the more versatile tool 
p::d::CellDataTransfer and consider the following:

   - Prepare a data container like `Vector>` where the outer 
   layer represents each cell in the current mesh, and the inner layer 
   corresponds to the dof values inside each cell.
   - Prepare data for the updated grid on the old grid. 
  - On already activated cells, store dof values with 
  `cell->get_interpolated_dof_values()`.
  - On all other cells, store an empty container.
  - Register your data container for and execute coarsening and 
   refinement.
   - Interpolate the old solution on the new mesh.
   - Initialize your new solution vector with invalid values 
  `std::numeric_limits::infinity`.
  - On previously activated cells, extract the stored data with 
  `cell->set_dof_values_by_interpolation()`.
  - Skip all other cells which only have an empty container.
   - On non-ghosted solution vectors, call `compress(VectorOperation::min)` 
   to get correct values on ghost cells.

This leaves you with a correctly interpolated solution on the new mesh, 
where all newly activated dofs have the value `infinity`.

You can now (and not earlier!!!) assign the values you have in mind for the 
newly activated dofs. You may want to exchange data on ghost cells once 
more with `GridTools::exchange_cell_data_to_ghosts()`.

This is a fairly new feature in the library for a very specific use case, 
so there is no dedicated class for transferring solutions across finite 
element activation yet. If you successfully manage to make this work, would 
you be willing to turn this into a class for the deal.II library?

Marc
On Wednesday, December 30, 2020 at 8:22:59 AM UTC-7 k.d...@gmail.com wrote:

> Hi all,
> Thank you for your reply. 
> Let me explain what I am trying to do and why. 
> I want to solve a transient heat transfer problem of the additive 
> manufacturing (AM) process. In AM processes, metal powder is deposited in 
> layers, and then a laser source scans each layer and melts and bonds the 
> powder to the layer underneath it. To simulate this layer by layer process, 
> I want to start with a grid that covers all the layers, but initially, only 
> the bottom-most layer is active and all other layers are inactive, and then 
> as time progresses, I want to activate one layer at a time. I read on this 
> mailing list that cell "birth" or "activation" can be done by changing a 
> cell from FE_Nothing to FE_Q using p-refinement. I was trying to keep all 
> cells of the grid initially to FE_Nothing except the bottom-most layer. And 
> then convert one layer at a time to FE_Q. My questions are:
> 1. Does this make sense? 
> 2. I have to do two other things for this to work. (a) When a cell becomes 
> FE_Q from FE_Nothing, dofs that are activating for the 1st time, I need to 
> apply a non-zero initial value to those dofs. This is to simulation the 
> metal powder deposited at a specified temperature,. e.g. room temperature. 
> (b) the dofs, that were shared between a FE_Q and FE_Nothing cells before 
> the p-refinement and now shared between FE_Q and FE_Nothing cells after 
> refinement, should retrain their values from before the refinement. 
>
> Another way to simulation this process would be to use a cell "awaking" 
> process, instead of the cell "birth". I keep call cells FE_Q but apply a 
> very low diffusivity to the cells of the layers that are not yet "awake". 
> But this way, I have to solve a larger system in all time steps. I was 
> hoping to save some computation time, by only forming a system consists of 
> cells that are in the "active" layers only. 
>
> Please let me if this makes sense? Is there any other method in deal.ii 
> that can simulation such a process? 
> Thank you very much and happy holidays.
> Kaushik 
>
>
> On Tue, Dec 29, 2020 at 12:26 PM Wolfgang Bangerth  
> wrote:
>
>> On 12/28/20 5:11 PM, Marc Fehling wrote:
>> > 
>> > In case a FE_Nothing has been configured to dominate, the solution 
>> should be 
>> > continuous on the interface if I understood correctly, i.e., zero on 
>> the face. 
>> > I will write a few tests to see if this is actually automatically the 
>> case in 
>> > user applications. If so, this check for domination will help other 
>> users to 
>> > avoid this pitfall :)
>> > 
>>
>> More tests = more better :-)
>> Cheers
>>   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

[deal.II] deal.II Newsletter #147

['Rene Gassmoeller' via deal.II User Group]

Hello everyone!

This is deal.II newsletter #147.
It automatically reports recently merged features and discussions about the 
deal.II finite element library.


## Below you find a list of recently proposed or merged features:

#11449: Convert step-2 to a simplex test (proposed by peterrum) 
https://github.com/dealii/dealii/pull/11449

#11448: Added test for FENothing and FEQ constraints. (proposed by marcfehling) 
https://github.com/dealii/dealii/pull/11448

#11447: Utility grid generator function `hyper_line` for tests. (proposed by 
marcfehling) https://github.com/dealii/dealii/pull/11447

#11446: Update grid/cell_id_08 output. (proposed by marcfehling) 
https://github.com/dealii/dealii/pull/11446

#11445: Minor change to CellWeights exception handling. (proposed by 
marcfehling) https://github.com/dealii/dealii/pull/11445

#11444: Refactor the HDF5 abort code. (proposed by drwells) 
https://github.com/dealii/dealii/pull/11444

#11443: Make sure every particle always has a property pool associated with it. 
(proposed by bangerth) https://github.com/dealii/dealii/pull/11443

#11442: Make a test more verbose. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/11442

#11441: Avoid warning about unused variable in particle.h (proposed by 
masterleinad; merged) https://github.com/dealii/dealii/pull/11441

#11440: Avoid MPI calls in the HDF5 destructors when there is an exception - 
DataSet (proposed by dangars; merged) 
https://github.com/dealii/dealii/pull/11440

#11439: Fix warning due to boost (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11439

#11438: Fix boost include file. (proposed by luca-heltai; merged) 
https://github.com/dealii/dealii/pull/11438

#11437: Clean up FE_Nothing::get_name(). (proposed by drwells; merged) 
https://github.com/dealii/dealii/pull/11437

#11436: Add test for geometric global-coarsening multigrid (proposed by 
peterrum; merged) https://github.com/dealii/dealii/pull/11436

#11434: Rename file in test folder (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11434

#11433: Add test for p-multigrid (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11433

#11432: MatrixFreeTools::compute_diagonal() make class const (proposed by 
peterrum; merged) https://github.com/dealii/dealii/pull/11432

#11431: CMake: Fix Ginkgo linkage (proposed by tamiko; merged) 
https://github.com/dealii/dealii/pull/11431

#11430: Better document FE_Nothing. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/11430

#11429: Generalize MGTwoLevelTransfer for simplex meshes (proposed by peterrum) 
https://github.com/dealii/dealii/pull/11429

#11428: FETools::get_projection_matrix for simplices (proposed by peterrum) 
https://github.com/dealii/dealii/pull/11428

#11427: Fix tests with output files with same name (proposed by peterrum; 
merged) https://github.com/dealii/dealii/pull/11427

#11426: Add utility functions to ReferenceCell namespace (proposed by peterrum; 
merged) https://github.com/dealii/dealii/pull/11426

#11425: Fix output of grid/cell_id_08 (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11425

#11424: Introduce ReferenceCell::compute_orientation and 
::permute_according_orientation (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11424

#11423: Introduce ReferenceCell::unit_tangential_vectors() (proposed by 
peterrum; merged) https://github.com/dealii/dealii/pull/11423

#11422: Work on compute_embedding_matrices() for simplex (proposed by peterrum) 
https://github.com/dealii/dealii/pull/11422

#11421: Fix formatting (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11421

#11420: Test alternative path in QProjector (proposed by peterrum) 
https://github.com/dealii/dealii/pull/11420

#11419: Pyramid/wedge: diverse fixes (proposed by peterrum; merged) 
https://github.com/dealii/dealii/pull/11419

#11418: step-27: Updated to new interface. (proposed by marcfehling) 
https://github.com/dealii/dealii/pull/11418

#11417: CI: remove unnecessary OSX setup steps (proposed by tjhei; merged) 
https://github.com/dealii/dealii/pull/11417

#11416: Simplex/pyramid/wedge: add DG test (proposed by peterrum) 
https://github.com/dealii/dealii/pull/11416

#11415: Ignore deprecated-volatile warning when disabling exta diagnostic 
(proposed by Rombur) https://github.com/dealii/dealii/pull/11415

#11414: Raviart-Thomas orientation and permutation issue (proposed by konsim83) 
https://github.com/dealii/dealii/pull/11414

#11413: Ginkgo Updates - 2nd try (proposed by tamiko; merged) 
https://github.com/dealii/dealii/pull/11413

#11412: Silence unused parameter warning when building with C++17 (proposed by 
masterleinad; merged) https://github.com/dealii/dealii/pull/11412

#11411: make constructors of Table and AlignedVector explicit (proposed by 
tjhei; merged) https://github.com/dealii/dealii/pull/11411

#11410: Revert "Merge pull request #11402 from pr

Re: [deal.II] Re: Modifying shape function data in MatrixFree

[Jean-Paul Pelteret]

Hi David,

One last point I don't understand right now is the magnitude of the 
machine precision you talk about:


> You'll see that they all exhibit quadratic convergence, with only a 
slight difference in the convergence history that can probably be 
attributed to accumulated round-off errors. Nevertheless, for all 
implementations the absolute residual at the end of each timestep can 
be pushed down to near machine precision with only a few Newton 
iterations.
I was just indicating that if you take sufficient Newton steps then the 
out-of-balance forces are of the order of machine precision, and its 
unlikely that the solution could be improved by any appreciable margin. 
Here's what happens when 4 Newton steps are used in the first timestep 
with the "tau" parameterisation:


Timestep 1 @ 0.1s
*** Using tau=tau(b) parameterisation ***
...
  4  ---  ASM_SYS  CONVERGED!
...
Absolute errors:
Displacement:    1.447e-09
Force:         1.018e-09

In this case, we've hit the early-exit criteria for the Newton scheme 
after 4 iterations. But for the "P" parameterisation, we miss this 
criterion by a small margin and take one more step.


Timestep 1 @ 0.1s
*** Using P=P(F) parameterisation ***
...
  5  ---  ASM_SYS  CONVERGED!
...
Absolute errors:
Displacement:    1.390e-12
Force:         8.075e-15

After that, the force residual (i.e. the norm of the "u" component of 
the RHS vector) is much smaller. But the actual quality of the solution, 
or the resulting stress field, probably haven't been improved by any 
significant margin (at least, not anything that impacts the scenario 
being simulated).


The absolute errors at the end of each iteration are of the order 
\mathcal{O}(10^{-10}). Also, if I compare my two-point residual 
assembly using matrix-free with the 'usual' spatial residual assembly 
and compute the l2 norm of the difference, i.e. l2_norm(r_mf - 
r-_tau), I obtain a difference of the order \mathcal{O}(10^{-10}). But 
from what I learnt, double precision should be at least accurate up to 
~10^{-15}. There are five orders of magnitude in between. I probably 
miss something here in between. Any idea? 


I understand that your description implies that you're comparing the 
residuals for a timestep where both schemes have taken the same number 
of Newton iterations. But nevertheless, I'm not quite sure that there's 
any value in comparing these two quantities. They are both discrete 
residuals for the approximate solution, and should both converge to zero 
norm. So neither one is actually the right reference value -- zero is. 
I'd explain this difference by the fact that both schemes use totally 
different code paths, so this is not an apples-to-apples comparison.  
The evolution of the solution for the MF and MB solvers might be such 
that the DoFs where the (very small) residuals remain at the end of each 
timestep differ. I'm not sure that one can say that this is actually 
means anything -- does it really matter that the displacement at one 
vertex differs by minuscule amount? Probably not. But even a small shift 
in the solution changes "where" the error in the residual resides, and 
your comparison is measuring that.


Best,
Jean-Paul

On 31.12.20 12:09, 'David' via deal.II User Group wrote:

Hi Jean-Paul,


> The indexing is correct - I've taken that part of the code from some 
other codes that I've fully verified. But I thought it would still be 
a useful exercise to re-implement the assembly routine for the other 
two parameterisations. Maybe this would further convince you of its 
correctness.


Convinced ;)

> I've attached a quickly modified version of step-44 that does this 
(see the functions assemble_system_one_cell_tau(), 
assemble_system_one_cell_S() and assemble_system_one_cell_P()), and 
uses the aforementioned push/pull operations to transform the stresses 
and material tangents. (I didn't feel like reimplementing the 
constitutive laws for the different parameterisations -- I leave that 
as an exercise to you, if you feel so inclined.) I've also attached 
some result logs for the default configuration for step-44 plus one 
additional global refinement step.


ah very nice, thanks again.

> That's great! To me, this implies that you fixed a bug. Ultimately 
the RHS that is computed via any of these parameterisations should, in 
principle, lead to exactly the same result. They are equivalent, and 
it is sometimes simply convenience (in terms of material laws, 
numerical efficiency,...) that dictates which you might choose. They 
are all expressing the weak form of the balance of linear momentum 
(i.e. the different power conjugate pairs collectively quantify the 
same mechanical power), and can happily be transformed from one to the 
other by exploiting the relationships between the various stresses and 
strains. This implies that the exact linearisation, even if its 
expressed in terms of some other stress-strain measures, will in fact 
be the linearisation of any one of these thr

Re: [deal.II] Modifying shape function data in MatrixFree

[luca.heltai]

> On 31 Dec 2020, at 12:09, 'David' via deal.II User Group 
>  wrote:
> 
> 
> The absolute errors at the end of each iteration are of the order 
> \mathcal{O}(10^{-10}). Also, if I compare my two-point residual assembly 
> using matrix-free with the 'usual' spatial residual assembly and compute the 
> l2 norm of the difference, i.e. l2_norm(r_mf - r-_tau), I obtain a difference 
> of the order \mathcal{O}(10^{-10}). But from what I learnt, double precision 
> should be at least accurate up to ~10^{-15}. There are five orders of 
> magnitude in between. I probably miss something here in between. Any idea? 

Well, you have to keep in mind that there are two more ingredients to keep into 
account: 

1. the Jacobian condition number (easily in the order of 10^3 — 10^6)
2. the number of operations that lead to the residuals (i.e., how many sums 
your are eventually doing)

since errors accumulate, 10^{-10} is actually close to machine precision if you 
have around 10^4/10^5 dofs. 

Moreover, you are computing the l2_norm, i.e., sum over all dofs  of the 
*square* of the difference, and then you are taking the square root. 

If you make 1e-15 error et each step in the square difference, propagate that 
for ndofs sums, simply by taking the square root of the resulting error, you’d 
get an expected error close to 1^{-8}, so you’re results are essentially 
machine precision.

L.

-- 
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/9D1619AB-95F5-4BE8-9360-162CBB69D9C5%40gmail.com.

Re: [deal.II] Re: Modifying shape function data in MatrixFree

['David' via deal.II User Group]

Hi Jean-Paul,


> The indexing is correct - I've taken that part of the code from some 
other codes that I've fully verified. But I thought it would still be a 
useful exercise to re-implement the assembly routine for the other two 
parameterisations. Maybe this would further convince you of its 
correctness. 

Convinced ;)

> I've attached a quickly modified version of step-44 that does this (see 
the functions assemble_system_one_cell_tau(), assemble_system_one_cell_S() 
and assemble_system_one_cell_P()), and uses the aforementioned push/pull 
operations to transform the stresses and material tangents. (I didn't feel 
like reimplementing the constitutive laws for the different 
parameterisations -- I leave that as an exercise to you, if you feel so 
inclined.) I've also attached some result logs for the default 
configuration for step-44 plus one additional global refinement step.

ah very nice, thanks again.

> That's great! To me, this implies that you fixed a bug. Ultimately the 
RHS that is computed via any of these parameterisations should, in 
principle, lead to exactly the same result. They are equivalent, and it is 
sometimes simply convenience (in terms of material laws, numerical 
efficiency,...) that dictates which you might choose. They are all 
expressing the weak form of the balance of linear momentum (i.e. the 
different power conjugate pairs collectively quantify the same mechanical 
power), and can happily be transformed from one to the other by exploiting 
the relationships between the various stresses and strains. This implies 
that the exact linearisation, even if its expressed in terms of some other 
stress-strain measures, will in fact be the linearisation of any one of 
these three expressions of the residual associated with displacement DoFs. 
The linearisation can be similarly computed for one parameterisation and 
transformed into the others. Its a nice exercise to go though, if you have 
the time and patience to do so.

Ah ok, I was not aware of that! In your first answer, you wrote "So you may 
as well start off by parameterising the problem in terms of F and P, and 
naturally you have to adjust the linearisation as well." I thought I need 
to adjust the linearization in any case, when I start to rephrase the 
residual assembly, but I guess you talked about an entirely rephrased 
problem (i.e., residual and linearization), which requires an adaption of 
the linearization. I checked now the convergence behavior of my two-point 
residual assembly in combination with the spatial linearization and I 
observe indeed quadratic convergence (though there are some deviations 
depending on the particular iteration). Rephrasing the linearization is 
therefore as an exercise left for the next pandemic.
One last point I don't understand right now is the magnitude of the machine 
precision you talk about:

> You'll see that they all exhibit quadratic convergence, with only a 
slight difference in the convergence history that can probably be 
attributed to accumulated round-off errors. Nevertheless, for all 
implementations the absolute residual at the end of each timestep can be 
pushed down to near machine precision with only a few Newton iterations.

The absolute errors at the end of each iteration are of the order \mathcal
{O}(10^{-10}). Also, if I compare my two-point residual assembly using 
matrix-free with the 'usual' spatial residual assembly and compute the l2 
norm of the difference, i.e. l2_norm(r_mf - r-_tau), I obtain a difference 
of the order \mathcal{O}(10^{-10}). But from what I learnt, double 
precision should be at least accurate up to ~10^{-15}. There are five 
orders of magnitude in between. I probably miss something here in between. 
Any idea? 

Regards,
David

Jean-Paul Pelteret schrieb am Mittwoch, 30. Dezember 2020 um 21:27:30 UTC+1:

> HI David,
>
> You're welcome! I'm glad that you found the explanation useful :-)
>
>
> the required changes for the tangent operator seem acceptable (assuming 
> the indexing in the get_HH() function does the right thing, otherwise it's 
> probably not acceptable)
>
>
> The indexing is correct - I've taken that part of the code from some other 
> codes that I've fully verified. But I thought it would still be a useful 
> exercise to re-implement the assembly routine for the other two 
> parameterisations. Maybe this would further convince you of its 
> correctness. 
>
> I've attached a quickly modified version of step-44 that does this (see 
> the functions assemble_system_one_cell_tau(), assemble_system_one_cell_S() 
> and assemble_system_one_cell_P()), and uses the aforementioned push/pull 
> operations to transform the stresses and material tangents. (I didn't feel 
> like reimplementing the constitutive laws for the different 
> parameterisations -- I leave that as an exercise to you, if you feel so 
> inclined.) I've also attached some result logs for the default 
> configuration for step-44 plus one additional global refinement