[deal.II] Re: online deal.II workshop: Friday, June 18

['David' via deal.II User Group]

Hi Timo,

short question: will the workshop afterwards be available on YouTube 
analogous to the last workshop?

Regards, 
David
On Wednesday, 2 June 2021 at 17:36:08 UTC+2 Timo Heister wrote:

> Hi all,
>
> We would like to announce a one-day deal.II workshop on June 18, 2021
> with several talks about recent developments in the library, the 9.3
> release, and interesting user projects. The talks will be available
> live using zoom and available as a video at a later time.
>
> For more information including a schedule, please see:
> https://dealii.org/workshop-2021/
>
> Information on how to join will be made available closer to the date.
> We hope to see you there!
>
> Best,
> Timo and the deal.II developers
>
> -- 
> 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/084cee70-63c3-43ed-9202-574f8df34aa6n%40googlegroups.com.

Re: [deal.II] Unique FE evaluation on arbitrary points

['David' via deal.II User Group]

Hi Peter,

thanks for your quick reply!

> A quick (and dirty) solution would be to evaluate and perform the 
communication non-uniquely and pick out one of the received values at the 
end. See also VectorTools::evaluate_at_points() and the flag 
VectorTools::EvaluationFlags::insert: although intentionally not documented 
this function will take the first value of the process with the lowest rank.

I have not yet looked in such a detail at the functions, but if I have the 
solution values again on multiple processes I need to decide again on each 
process. However, let me first take a closer look at the functions, I just 
realized that there is also a `ConsensusAlgorithm` namespace which might 
offer some remedies. 

> Anyway, I am happy to discuss how to implement the needed feature. The 
code that would need to be touched is: 
https://github.com/dealii/dealii/blob/654b14481eeff7e744297025ea2078ab115611ac/source/grid/grid_tools.cc#L5905-L5909,
 
i.e., a consistent modification of send_components and recv_components. I 
am happy to write a draft if you want.

All in all, I don't need a mature solution right now and having something 
more compact inside the library would be of course highly appreciated. I 
could contribute something in this direction, if that's what you were 
thinking of. Let me add that it will probably take more time than you 
think. Looking forward to the draft.

David
peterr...@gmail.com schrieb am Mittwoch, 12. Mai 2021 um 12:30:21 UTC+2:

> Hi David,
>
> > I want to evaluate each point out of the cloud only once in case of a 
> distributed triangulation.
>
> A quick (and dirty) solution would be to evaluate and perform the 
> communication non-uniquely and pick out one of the received values at the 
> end. See also VectorTools::evaluate_at_points() and the flag 
> VectorTools::EvaluationFlags::insert: although intentionally not documented 
> this function will take the first value of the process with the lowest rank.
>
> > I would probably need to go a step further and implement a consensus 
> algorithm in order to decide only for one specific process, which seems not 
> trivial to me.
>
> You are right. One would need an additional communication step (with or 
> without CA) to make the evaluation/communication more efficient. I was 
> thinking about introducing an additional flag to the constructor of 
> RemotePointEvaluation, but have not implemented it yet, since in our use 
> cases we either need the average values (DG) or we set up the communication 
> pattern in every time step due to some very dynamic systems (multi phase 
> flow).
>
> Anyway, I am happy to discuss how to implement the needed feature. The 
> code that would need to be touched is: 
> https://github.com/dealii/dealii/blob/654b14481eeff7e744297025ea2078ab115611ac/source/grid/grid_tools.cc#L5905-L5909,
>  
> i.e., a consistent modification of send_components and recv_components. I 
> am happy to write a draft if you want.
>
> Hope this helps,
> Peter
> On 12.05.21 09:08, 'David' via deal.II User Group wrote:
>
> Dear all, 
>
> I'm currently looking for a way to evaluate my global solution on an 
> arbitrary cloud of points and stumbled across the (new) class 
> `Utilities::MPI::RemotePointEvaluation`. In general, the class does 
> exactly what I want with one particular exception; I want to evaluate each 
> point out of the cloud only once in case of a distributed triangulation. If 
> I understand it correctly this property is fulfilled in case the boolean 
> function `is_map_unique()` returns true, but I don’t see any option to 
> filter or tell the class that I actually want the map to be unique.
>
> Assuming that the class cannot handle this case I tried to pre-filter the 
> vertices by hand: For me, the most natural choice to fulfill the unique 
> evaluation of the points is to assign them to a process, so that I only 
> work on ‘locally owned vertices’. I compute a bounding box using the 
> iterator filter `LocallyOwnedCell()`, filter the received cloud of points 
> on the process and finally plug them into the `RemotePointEvaluation`. 
> However, this does not work in case points of the cloud are located at the 
> boundary of the bounding box, because they will technically be part of both 
> bounding boxes of the involved processes. I would probably need to go a 
> step further and implement a consensus algorithm in order to decide only 
> for one specific process, which seems not trivial to me. 
>
> Before starting now to implement a more advanced solution I would like to 
> ask if there is an obvious way I’m not aware of. Is it possible to enforce 
> a unique evaluation using one of the ‘PointEvaluation’ classes or is it 
> possible to let deal.II decide for an owner of an arbitrary vertex? I have 
> also seen

[deal.II] Unique FE evaluation on arbitrary points

['David' via deal.II User Group]

Dear all, 

I'm currently looking for a way to evaluate my global solution on an 
arbitrary cloud of points and stumbled across the (new) class 
`Utilities::MPI::RemotePointEvaluation`. In general, the class does 
exactly what I want with one particular exception; I want to evaluate each 
point out of the cloud only once in case of a distributed triangulation. If 
I understand it correctly this property is fulfilled in case the boolean 
function `is_map_unique()` returns true, but I don’t see any option to 
filter or tell the class that I actually want the map to be unique.

Assuming that the class cannot handle this case I tried to pre-filter the 
vertices by hand: For me, the most natural choice to fulfill the unique 
evaluation of the points is to assign them to a process, so that I only 
work on ‘locally owned vertices’. I compute a bounding box using the 
iterator filter `LocallyOwnedCell()`, filter the received cloud of points 
on the process and finally plug them into the `RemotePointEvaluation`. 
However, this does not work in case points of the cloud are located at the 
boundary of the bounding box, because they will technically be part of both 
bounding boxes of the involved processes. I would probably need to go a 
step further and implement a consensus algorithm in order to decide only 
for one specific process, which seems not trivial to me. 

Before starting now to implement a more advanced solution I would like to 
ask if there is an obvious way I’m not aware of. Is it possible to enforce 
a unique evaluation using one of the ‘PointEvaluation’ classes or is it 
possible to let deal.II decide for an owner of an arbitrary vertex? I have 
also seen the function ‘find_active_cell_around_point()’ which might help, 
but it might end up in similar problem as well.

Thanks in advance and 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c29fcdf5-59bf-450f-9de6-d9f9bddd3c68n%40googlegroups.com.

Re: [deal.II] Visualizing higher-order cells in 3D (step-67)

['David' via deal.II User Group]

Hi Peter,

>David, could you try out the following PR:
https://github.com/dealii/dealii/pull/11784

I tried it and the PR doesn't fix it for me.

> Not sure if the issue there is related. If not, could you post a vtk/vtu
file and I'll read it with my Paraview version.

Very good idea. I attached an affected file.
Can you reproduce the issue with your paraView version?

Regards,
David

Am Di., 23. Feb. 2021 um 16:08 Uhr schrieb peter rum :

> David, could you try out the following PR:
> https://github.com/dealii/dealii/pull/11784
> Not sure if the issue there is related. If not, could you post a vtk/vtu
> file and I'll read it with my Paraview version.
>
> Peter
>
> On Tuesday, 23 February 2021 at 09:07:57 UTC+1 daschn...@googlemail.com
> wrote:
>
>> Hi Wolfgang,
>>
>> > I have to admit that I don't know how to address this because I don't
>> know nor
>> > use Paraview. As mentioned before, it would be useful to obtain the
>> smallest
>> > possible testcase that reproduces the problem -- minimal number of
>> processors,
>> > minimal number of cells.
>>
>> I already tried to answer this question above since Alexander asked the
>> same question (not sure why his message has been deleted):
>>
>> > Yes, I can reproduce it on 2 procs and I can also simplify the mesh to
>> a size of four cells on a rectangle. The problem vanishes if I select
>> linear shape functions or zero nonlinear subdivisions in paraView (which is
>> probably the same).
>>
>> [image: coarse_mesh_distorted.png]
>>
>> > If I increase the nonlinear subdivision setting in paraVIew to four (as
>> opposed to the default one), it looks as follows:
>>
>> [image: coarse_mesh_four_subd.png]
>>
>>
>> Does it help?
>> Regards,
>> David
>>
>> Wolfgang Bangerth schrieb am Montag, 22. Februar 2021 um 20:10:05 UTC+1:
>>
>>> On 2/18/21 6:39 AM, 'David' via deal.II User Group wrote:
>>> >
>>> > the problem here is independent of any data set. So, you can also see
>>> it in
>>> > paraView in the 'solid color' block. It seems to be an artificial
>>> curvature
>>> > within the element.
>>> > To my best knowledge it is not required to update any geometry related
>>> data.
>>>
>>> I have to admit that I don't know how to address this because I don't
>>> know nor
>>> use Paraview. As mentioned before, it would be useful to obtain the
>>> smallest
>>> possible testcase that reproduces the problem -- minimal number of
>>> processors,
>>> minimal number of cells.
>>>
>>> Best
>>> Wolfgang
>>>
>>> --
>>> 
>>> 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 a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/O059tF8sy3o/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/1d13c760-0972-450f-a553-b26f2b2fb7den%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/1d13c760-0972-450f-a553-b26f2b2fb7den%40googlegroups.com?utm_medium=email_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/CAK3rKvE0q--moxsAUHj2GkLHUpSyfsjniTVvW3M01peKZ6%2BfYQ%40mail.gmail.com.


solution_000.vtu
Description: Binary data

Re: [deal.II] Visualizing higher-order cells in 3D (step-67)

['David' via deal.II User Group]

Hi Daniel,

the problem here is independent of any data set. So, you can also see it in 
paraView in the 'solid color' block. It seems to be an artificial curvature 
within the element. 
To my best knowledge it is not required to update any geometry related data.
What do you think?

Regards,
David

d.arnd...@gmail.com schrieb am Freitag, 12. Februar 2021 um 23:09:59 UTC+1:

> David,
>
> For me the most likely suspect without looking into the code at all is a 
> missing update_ghost_values.
>
> Best,
> Daniel
>
> Am Di., 9. Feb. 2021 um 17:04 Uhr schrieb Alexander :
>
>> David
>> i believe that in order to proceed, one will have to simplify this 
>> further.
>>
>> What is the minimum number of processes this happens? Can you reproduce 
>> it for 2 procs? 
>>
>> Additionally, can the mesh size be reduced to a minimum (ideally a 
>> handful of cells, say 4-8)?
>>
>> Alexander
>>
>> On Monday, February 8, 2021 at 10:20:35 AM UTC+1 daschn...@googlemail.com 
>> wrote:
>>
> Hi Wolfgang,
>>>
>>> > David -- can you be more explicit about what "a bit odd" means in 
>>> this case? 
>>> > What are we looking at, and why do you think this is wrong? 
>>>
>>> Certainly. We are looking at solution_000.vtu (zx-plane) of step-67, 
>>> where I changed the testcase to 1 (cylinder in channel) and the 
>>> dimension to 3D, 
>>> <https://github.com/dealii/dealii/blob/master/examples/step-67/step-67.cc#L77-L78>
>>>  
>>> computed on 4 cores.
>>>
>>> What I mean by 'a bit odd' are in particular the wrinkles on the 
>>> surface:
>>> [image: close-up.png]
>>>
>>>
>>> Running the same case in serial looks as follows.
>>> [image: serial.png]
>>>
>>> The surface is smooth as usual. I would expect the visualization of the 
>>> serial and the parallel case to be the same.
>>>
>>> Best regards,
>>> David
>>> Wolfgang Bangerth schrieb am Montag, 8. Februar 2021 um 05:12:42 UTC+1:
>>>
>>>> On 2/7/21 2:12 AM, 'David' via deal.II User Group wrote: 
>>>> > , 
>>>> > 
>>>> > I'm running a 3D case using the 'write-higher-order-cells' flag and 
>>>> the 
>>>> > 'write-vtu-in-parallel' function. The output writing is quite similar 
>>>> to the 
>>>> > way step-67 handles it. However, the output of my 3D data sets looks 
>>>> a bit odd 
>>>> > when running the case in parallel. 
>>>>
>>>> David -- can you be more explicit about what "a bit odd" means in this 
>>>> case? 
>>>> What are we looking at, and why do you think this is wrong? 
>>>>
>>>> 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+un...@googlegroups.com.
>
>
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/aaf0032a-63d8-4709-b913-d0932809f3b3n%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/aaf0032a-63d8-4709-b913-d0932809f3b3n%40googlegroups.com?utm_medium=email_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/0ed9d893-294b-405b-9ab9-a712438e657cn%40googlegroups.com.

Re: [deal.II] Understanding loops in matrix-free

['David' via deal.II User Group]

..just noted that the order is exactly the point in the loop framework.

I cannot assume a certain order in my boundary data (or I need to choose an 
order corresponding to the loop function).

Special thanks to @peterrum for pointing this out. This thread is therefore 
resolved.

Thank you very much!

David schrieb am Donnerstag, 18. Februar 2021 um 14:35:52 UTC+1:

> Hi Martin,
>
> thank you for your reply.
> before going into detail and providing more information let me ask again 
> in order to avoid a misunderstanding. You say
>
> > apart from when exactly the cells and faces are scheduled, i.e., 
> possibly a change in the order of roundoff errors.
>
> As I said, my boundary condition is in an ordered vector (according to the 
> face loop) and in my 'local_apply_boundary_face' I explicitly iterate over 
> the boundary data. In pseudo code, it would look like this:
>
> int iterator = 0
> for (boundary_face_loop)
> {
>  // apply boundary condition contained in
>  my_ordered_boundary_data[ iterator ]
> ++ iterator;
> }
>
> If I change now the order (when faces are scheduled), I would of course 
> apply different data to different faces. I'm not sure if I exactly 
> understand the meaning of 'when exactly cells and faces are scheduled'. 
> But I assume the implementation/ordering is not a problem.
> I have another idea how to check it and will come back to this topic after 
> investigating.
>
> Regards,
> David
>
> Martin Kronbichler schrieb am Samstag, 13. Februar 2021 um 15:53:11 UTC+1:
>
>> Dear David,
>>
>> Indeed the two approaches should give the same answer; apart from when 
>> exactly the cells and faces are scheduled, i.e., possibly a change in the 
>> order of roundoff errors. Thus, I suspect there is something strange going 
>> on and could be a bug either in deal.II or in how you set up some 
>> particular ingredients. Can you share the code or explain a bit what 
>> exactly is used? We have continuous elements, but which degree? Is it a 
>> system? How many cells? Are you using threads, i.e., is any of the task 
>> parallel settings turned on? From what you say I believe you are running in 
>> serial, right? So it seems a pretty small problem, which we should be able 
>> to investigate rapidly.
>>
>> One thing you could do is to look into the "face_range" that you obtain 
>> when the algorithm calls back into local_apply_boundary_face and compare 
>> that with the range that you manually construct in your first version? I 
>> wonder if there are some parts of the loop we are missing or running twice.
>>
>> Best,
>> Martin
>> On 06.02.21 19:18, 'David' via deal.II User Group wrote:
>>
>>
>> Sorry for messing up the topic. I should be Understanding loops in 
>> matrix-free. I wanted to insert a figure of the source code rather than the 
>> google groups formatting and it didn't work for some reason.
>> David schrieb am Samstag, 6. Februar 2021 um 19:15:07 UTC+1:
>>
>>> Hi there, 
>>>
>>> I'm currently trying to pack my cell operations into one of the 
>>> matrix-free loop functions.
>>> In my first version, I implemented the loops manually by using (sorry 
>>> for the odd formatting):
>>> ```
>>> local_apply_cell(*matrix_free,
>>>  system_rhs,
>>>  source_vector,
>>>  std::make_pair(0,
>>> 
>>> matrix_free->n_cell_batches()));
>>> local_apply_boundary_face(
>>>   *matrix_free,
>>>   system_rhs,
>>>   source_vector,
>>>   std::make_pair(mf_data_reference->n_inner_face_batches(),
>>>  mf_data_reference->n_inner_face_batches() +
>>>
>>> mf_data_reference->n_boundary_face_batches()));
>>> ```
>>> which works as desired. I don't have any operations on internal faces 
>>> and ghost exchange/compression doesn't matter for the rest of the question.
>>>
>>> In a second step, I replaced the manually implemented loops by the 
>>> matrix-free loop in the following way:
>>> ```
>>> matrix_free->template loop(
>>>   [&](const auto ,
>>>   auto &  dst,
>>>   const auto ,
>>>   const auto _range) {
>>> local_apply_cell(data, dst, src, cell_range);
>>>  

Re: [deal.II] Understanding loops in matrix-free

['David' via deal.II User Group]

Hi Martin,

thank you for your reply.
before going into detail and providing more information let me ask again in 
order to avoid a misunderstanding. You say

> apart from when exactly the cells and faces are scheduled, i.e., possibly 
a change in the order of roundoff errors.

As I said, my boundary condition is in an ordered vector (according to the 
face loop) and in my 'local_apply_boundary_face' I explicitly iterate over 
the boundary data. In pseudo code, it would look like this:

int iterator = 0
for (boundary_face_loop)
{
 // apply boundary condition contained in
 my_ordered_boundary_data[ iterator ]
++ iterator;
}

If I change now the order (when faces are scheduled), I would of course 
apply different data to different faces. I'm not sure if I exactly 
understand the meaning of 'when exactly cells and faces are scheduled'. 
But I assume the implementation/ordering is not a problem.
I have another idea how to check it and will come back to this topic after 
investigating.

Regards,
David

Martin Kronbichler schrieb am Samstag, 13. Februar 2021 um 15:53:11 UTC+1:

> Dear David,
>
> Indeed the two approaches should give the same answer; apart from when 
> exactly the cells and faces are scheduled, i.e., possibly a change in the 
> order of roundoff errors. Thus, I suspect there is something strange going 
> on and could be a bug either in deal.II or in how you set up some 
> particular ingredients. Can you share the code or explain a bit what 
> exactly is used? We have continuous elements, but which degree? Is it a 
> system? How many cells? Are you using threads, i.e., is any of the task 
> parallel settings turned on? From what you say I believe you are running in 
> serial, right? So it seems a pretty small problem, which we should be able 
> to investigate rapidly.
>
> One thing you could do is to look into the "face_range" that you obtain 
> when the algorithm calls back into local_apply_boundary_face and compare 
> that with the range that you manually construct in your first version? I 
> wonder if there are some parts of the loop we are missing or running twice.
>
> Best,
> Martin
> On 06.02.21 19:18, 'David' via deal.II User Group wrote:
>
>
> Sorry for messing up the topic. I should be Understanding loops in 
> matrix-free. I wanted to insert a figure of the source code rather than the 
> google groups formatting and it didn't work for some reason.
> David schrieb am Samstag, 6. Februar 2021 um 19:15:07 UTC+1:
>
>> Hi there, 
>>
>> I'm currently trying to pack my cell operations into one of the 
>> matrix-free loop functions.
>> In my first version, I implemented the loops manually by using (sorry for 
>> the odd formatting):
>> ```
>> local_apply_cell(*matrix_free,
>>  system_rhs,
>>  source_vector,
>>  std::make_pair(0,
>> 
>> matrix_free->n_cell_batches()));
>> local_apply_boundary_face(
>>   *matrix_free,
>>   system_rhs,
>>   source_vector,
>>   std::make_pair(mf_data_reference->n_inner_face_batches(),
>>  mf_data_reference->n_inner_face_batches() +
>>
>> mf_data_reference->n_boundary_face_batches()));
>> ```
>> which works as desired. I don't have any operations on internal faces and 
>> ghost exchange/compression doesn't matter for the rest of the question.
>>
>> In a second step, I replaced the manually implemented loops by the 
>> matrix-free loop in the following way:
>> ```
>> matrix_free->template loop(
>>   [&](const auto ,
>>   auto &  dst,
>>   const auto ,
>>   const auto _range) {
>> local_apply_cell(data, dst, src, cell_range);
>>   },
>>   [](const auto &, auto &, const auto &, const auto &) {},
>>   [&](const auto ,
>>   auto &  dst,
>>   const auto ,
>>   const auto _range) {
>> local_apply_boundary_face(data, dst, src, face_range);
>>   },
>>   system_rhs,
>>   source_vector,
>>   true,
>>   MatrixFree> VectorizedArrayType>::DataAccessOnFaces::
>> none,
>>   MatrixFree> VectorizedArrayType>::DataAccessOnFaces::
>> none);
>> ```
>>
>> However, the system

[deal.II] Understanding loops in matrix-free

['David' via deal.II User Group]

Sorry for messing up the topic. I should be Understanding loops in 
matrix-free. I wanted to insert a figure of the source code rather than the 
google groups formatting and it didn't work for some reason.
David schrieb am Samstag, 6. Februar 2021 um 19:15:07 UTC+1:

> Hi there,
>
> I'm currently trying to pack my cell operations into one of the 
> matrix-free loop functions.
> In my first version, I implemented the loops manually by using (sorry for 
> the odd formatting):
> ```
> local_apply_cell(*matrix_free,
>  system_rhs,
>  source_vector,
>  std::make_pair(0,
> 
> matrix_free->n_cell_batches()));
> local_apply_boundary_face(
>   *matrix_free,
>   system_rhs,
>   source_vector,
>   std::make_pair(mf_data_reference->n_inner_face_batches(),
>  mf_data_reference->n_inner_face_batches() +
>
> mf_data_reference->n_boundary_face_batches()));
> ```
> which works as desired. I don't have any operations on internal faces and 
> ghost exchange/compression doesn't matter for the rest of the question.
>
> In a second step, I replaced the manually implemented loops by the 
> matrix-free loop in the following way:
> ```
> matrix_free->template loop(
>   [&](const auto ,
>   auto &  dst,
>   const auto ,
>   const auto _range) {
> local_apply_cell(data, dst, src, cell_range);
>   },
>   [](const auto &, auto &, const auto &, const auto &) {},
>   [&](const auto ,
>   auto &  dst,
>   const auto ,
>   const auto _range) {
> local_apply_boundary_face(data, dst, src, face_range);
>   },
>   system_rhs,
>   source_vector,
>   true,
>   MatrixFree VectorizedArrayType>::DataAccessOnFaces::
> none,
>   MatrixFree VectorizedArrayType>::DataAccessOnFaces::
> none);
> ```
>
> However, the system starts diverging after a time (as opposed to the 
> 'manual' loop), i.e., the result is different.
>
> I should probably add as a comment that I use sorted boundary condition 
> according to the boundary cell batches so that I apply the first set of 
> data to the first boundary face batch and the second one to the second... 
> Though, I checked the ordering manually and the loop-function seems to work 
> in the same order on the batches as the manual loop.
>
> I also tried to set the vector zeroing on the destination vector to 
> 'false' (I use continuous elements) but it doesn't make any difference. 
> According to my understanding of the loop function, both implementations 
> should be equivalent. I probably miss something here.
>
> Maybe anyone else has an idea about the differences.
> Thanks in advance and 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/620dc4c7-4adf-47f7-8f43-cc9a3dd3e22en%40googlegroups.com.

[deal.II] Understanding

['David' via deal.II User Group]

Hi there,

I'm currently trying to pack my cell operations into one of the matrix-free 
loop functions.
In my first version, I implemented the loops manually by using (sorry for 
the odd formatting):
```
local_apply_cell(*matrix_free,
 system_rhs,
 source_vector,
 std::make_pair(0,
matrix_free->n_cell_batches()));
local_apply_boundary_face(
  *matrix_free,
  system_rhs,
  source_vector,
  std::make_pair(mf_data_reference->n_inner_face_batches(),
 mf_data_reference->n_inner_face_batches() +
   
mf_data_reference->n_boundary_face_batches()));
```
which works as desired. I don't have any operations on internal faces and 
ghost exchange/compression doesn't matter for the rest of the question.

In a second step, I replaced the manually implemented loops by the 
matrix-free loop in the following way:
```
matrix_free->template loop(
  [&](const auto ,
  auto &  dst,
  const auto ,
  const auto _range) {
local_apply_cell(data, dst, src, cell_range);
  },
  [](const auto &, auto &, const auto &, const auto &) {},
  [&](const auto ,
  auto &  dst,
  const auto ,
  const auto _range) {
local_apply_boundary_face(data, dst, src, face_range);
  },
  system_rhs,
  source_vector,
  true,
  MatrixFree::DataAccessOnFaces::
none,
  MatrixFree::DataAccessOnFaces::
none);
```

However, the system starts diverging after a time (as opposed to the 
'manual' loop), i.e., the result is different.

I should probably add as a comment that I use sorted boundary condition 
according to the boundary cell batches so that I apply the first set of 
data to the first boundary face batch and the second one to the second... 
Though, I checked the ordering manually and the loop-function seems to work 
in the same order on the batches as the manual loop.

I also tried to set the vector zeroing on the destination vector to 'false' 
(I use continuous elements) but it doesn't make any difference. According 
to my understanding of the loop function, both implementations should be 
equivalent. I probably miss something here.

Maybe anyone else has an idea about the differences.
Thanks in advance and 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/090ebbeb-4207-4342-accd-9a0ef111697fn%40googlegroups.com.

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

['David' via deal.II User Group]

 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 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: [deal.II] Re: Modifying shape function data in MatrixFree

['David' via deal.II User Group]

ed some result logs for the default 
> configuration for step-44 plus one additional global refinement step. 
> 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.
>
> Comparing the P-formulation to the rest, you can see that its actually 
> quite straight forward to implement, so hopefully converting it to matrix 
> free wouldn't pose too much of a challenge. Note that in both the P- and 
> S-formulations, all of the linearisation terms coupled to the displacement 
> also change slightly. The coupling terms derive from the definition of the 
> hydrostatic stress, which is expressed differently for each 
> parameterisation.
>
>
> As an interesting side note: If I remove my wrong 'symm_grad_Nx' 
> contribution from my shared code snippet, the remaining code corresponds 
> already to a rephrased formulation using the Kirchoff stress and 
> deformation gradient (as you have already noticed). Using this residual in 
> combination with the "tau-Jc" tangent from step-44 (one-field adopted) 
> leads to a convergence of the nonlinear scheme. So, it seems the 
> linearization is still "correct enough" for the modified (two-point) frame 
> of the residual. I will -of course- not rely on this half baked solution. 
>
>
> 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.
>
> Best,
> Jean-Paul
>
>
> On 30.12.20 15:07, 'David' via deal.II User Group wrote:
>
> Hi Jean-Paul,
>
> first of all, thank you very much for your comprehensive answer. It's a 
> little bit unfortunate that there is no way for the deformed test gradient 
> evaluation using matrix-free. However, your code snippets and the 
> explanation are traceable and 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). As always, I was a 
> little bit too ambitious with my time schedule for this task. I will 
> document your suggestions in an issue and try to implement the rephrased 
> problem in the medium run. In the end, the change is entirely performance 
> motivated, since the residual assembly (without matrx-free) is 
> comparatively expensive considering that it is only performed once per 
> nonlinear iteration. 
>
> As an interesting side note: If I remove my wrong 'symm_grad_Nx' 
> contribution from my shared code snippet, the remaining code corresponds 
> already to a rephrased formulation using the Kirchoff stress and 
> deformation gradient (as you have already noticed). Using this residual in 
> combination with the "tau-Jc" tangent from step-44 (one-field adopted) 
> leads to a convergence of the nonlinear scheme. So, it seems the 
> linearization is still "correct enough" for the modified (two-point) frame 
> of the residual. I will -of course- not rely on this half baked solution. 
>
> Best regards,
> David
>
> Jean-Paul Pelteret schrieb am Dienstag, 29. Dezember 2020 um 20:06:17 
> UTC+1:
>
>> Hi David,
>>
>> So as I'm sure that you know, if you want to assemble the linear system 
>> for quasi-static non-linear (hyper)elasticity with a referential DoFHandler 
>> (using no Eulerian mapping to transform the shape functions to the spatial 
>> setting) AND without doing a similar transformation by hand (as step-44 
>> does) then you are pretty much limited to using one of two 
>> parameterisations. E

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

['David' via deal.II User Group]

 case it makes a
> // difference) we transform it first.
> Tensor<4, dim> tmp1;
> 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)
> for (unsigned int K = 0; K < dim; ++K)
> tmp1[I][J][k][L] += get_F()[k][K] * H[I][J][K][L];
> Tensor<4, dim> HH_mixed;
> const Tensor<2, dim> I_ns = Physics::Elasticity::StandardTensors::I;
> 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)
> {
> // Push forward index 0 of material stress contribution
> for (unsigned int I = 0; I < dim; ++I)
> HH_mixed[i][J][k][L] += get_F()[i][I] * tmp1[I][J][k][L];
> // Add geometric stress contribution
> HH_mixed[i][J][k][L] += I_ns[i][k] * S[J][L];
> }
> return HH_mixed;
> }
> } 
> Naturally, the above could be simplified a bit for this fixed 
> parameterisation.
>
> With this, the (matrix-based) assembly would looks something like this for 
> the RHS
> // Excuse the messiness -- coding in an email client is not a good idea! 
> for (unsigned int q_point = 0; q_point < this->n_q_points; ++q_point) { const 
> Tensor<2,dim> P = lqph[q_point]->get_P(); 
> for (unsigned int i = 0; i < this->dofs_per_cell; ++i) {
> const unsigned int i_group = this->fe.system_to_base_index(i).first.first; 
> const Tensor<2,dim> dF_I = fe_values_ref[this->u_fe].gradient(i, q_point); 
>
> if (i_group == this->u_dof)
> data.cell_rhs(i) -= double_contract(dF_I, P) * JxW;
> ... // See the rest of step-44
>
> and the salient part of the tangent matrix assembly would be something like
> for (unsigned int q_point = 0; q_point < this->n_q_points; ++q_point) { // 
> Linearisation of P with respect to F; // contains both material and 
> geometric tangent contributions  const Tensor<4,dim> HH = 
> lqph[q_point]->get_HH(); 
>
> for (unsigned int i = 0; i < this->dofs_per_cell; ++i) {
> const unsigned int i_group = this->fe.system_to_base_index(i).first.first; 
> const Tensor<2,dim> dF_I = fe_values_ref[this->u_fe].gradient(i, q_point);
> for (unsigned int j = 0; j <= i; ++j)
> {
> const unsigned int j_group = this->fe.system_to_base_index(j).first.first;
> const Tensor<2,dim> dF_J = fe_values_ref[this->u_fe].gradient(j, q_point); 
>
> if ((i_group == j_group) && (i_group == this->u_dof))
> {
> cell_matrix(i, j) += scalar_product(dF_I, HH, dF_J) * JxW; } ... // See 
> the rest of step-44 
>
> So IIRC the call to phi.submit_gradient() that is bound to the 
> undeformed/non-mapped DoFHandler is effectively the same as testing against 
> what I've called dF_I in the above code. I think that the transformations 
> from tau->P and Jc->HH should make it simple enough to adjust what remains.
>
> I hope that this helps you implement what you're trying to do.
> Jean-Paul
>
>
>
> On 29.12.20 15:59, 'David' via deal.II User Group wrote:
>
> Maybe as an edit: what I currently do looks the following way:
>
> ```
> // Get gradient in reference frame
> const Tensor<2, dim, VectorizedArrayType> grad_u = phi.get_gradient(
> q_point); 
> // Compute deformation gradient
> const Tensor<2, dim, VectorizedArrayType> F = Physics::Elasticity::
> Kinematics::F(grad_u); 
> const SymmetricTensor<2, dim, VectorizedArrayType> b = Physics::Elasticity
> ::Kinematics::b(F); 
>
> const VectorizedArrayType det_F = determinant(F); 
> // Invert F
> const Tensor<2, dim, VectorizedArrayType> F_inv = invert(F); 
> const SymmetricTensor<2, dim, VectorizedArrayType> b_bar = Physics::
> Elasticity::Kinematics::b( Physics::Elasticity::Kinematics::F_iso(F)); 
> // Get tau from material description 
> SymmetricTensor<2, dim, VectorizedArrayType> tau; material->get_tau(tau, 
> det_F, b_bar, b);
>
> // Gradient itslef is included in the integrate call, apply push forward 
> with F^{-1}
> const Tensor<2, dim, VectorizedArrayType> symm_grad_Nx = symmetrize(F_inv
> );
> // Compute the result
> const Tensor<2, dim, VectorizedArrayType> result = symm_grad_Nx * Tensor<2
> , dim, VectorizedArrayType>(tau);
> // Apply an additional push forward with F^{-T}
> phi.submit_gradient(-result * transpose(F_inv), q_point);  // end loop 
> over quadrature points 
>
> // For each element
> phi.integrate(false, true); 
> ```
>
> It works, but I don't get the quadratic convergence order of the Newton 
> scheme anymore. Hence, I assume this is not sufficient.
> David schrieb am Dienstag, 2

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

['David' via deal.II User Group]

Maybe as an edit: what I currently do looks the following way:

```
// Get gradient in reference frame
const Tensor<2, dim, VectorizedArrayType> grad_u = phi.get_gradient(q_point
); 
// Compute deformation gradient
const Tensor<2, dim, VectorizedArrayType> F = Physics::Elasticity::
Kinematics::F(grad_u); 
const SymmetricTensor<2, dim, VectorizedArrayType> b = Physics::Elasticity::
Kinematics::b(F); 

const VectorizedArrayType det_F = determinant(F); 
// Invert F
const Tensor<2, dim, VectorizedArrayType> F_inv = invert(F); 
const SymmetricTensor<2, dim, VectorizedArrayType> b_bar = Physics::
Elasticity::Kinematics::b( Physics::Elasticity::Kinematics::F_iso(F)); 
// Get tau from material description 
SymmetricTensor<2, dim, VectorizedArrayType> tau; material->get_tau(tau, 
det_F, b_bar, b);

// Gradient itslef is included in the integrate call, apply push forward 
with F^{-1}
const Tensor<2, dim, VectorizedArrayType> symm_grad_Nx = symmetrize(F_inv);
// Compute the result
const Tensor<2, dim, VectorizedArrayType> result = symm_grad_Nx * Tensor<2, 
dim, VectorizedArrayType>(tau);
// Apply an additional push forward with F^{-T}
phi.submit_gradient(-result * transpose(F_inv), q_point);  // end loop over 
quadrature points 

// For each element
phi.integrate(false, true); 
```

It works, but I don't get the quadratic convergence order of the Newton 
scheme anymore. Hence, I assume this is not sufficient.
David schrieb am Dienstag, 29. Dezember 2020 um 12:41:35 UTC+1:

>
> Hi all,
>
> I'm currently trying to implement a vectorized variant of the residual 
> assembly as it is done in step-44/one of the corresponding code-gallery 
> examples using FEEvaluation objects in combination with matrix-free. 
>
> I was not able to find an appropriate solution for the given code line 
> and
>  
> the subsequent application for the assembly process: the major problem is 
> that the shape function gradients are defined with regard to the current 
> configuration, whereas my matrix-free object holds a mapping, which refers 
> to the initial configuration. The reference configuration mapping  of my 
> matrix-free object is desired as the final integration is actually 
> performed in the reference frame.
> A valid option is probably (I have not tried it) to use a matrix-free 
> object with a different mapping (MappingQEulerian) which directly evaluates 
> the shape gradients in the deformed configuration and use another 
> referential matrix-free object for the integration part. 
> The main reason to avoid this approach is that it would require a 
> reinitialization of the 'deformed' matrix-free object in each nonlinear 
> step and the mapping needs to be recomputed in each reinizialization. On 
> the other hand, the expensive reinitialization is not really required as 
> the mapping between the referential and the current configuration is known 
> due to the (inverse) deformation gradient. 
> Therefore, I'm looking for an opportunity to access the shape gradients 
> and perform a push forward similarly to the approach in step-44 in order to 
> evaluate the desired quantities. 
>
> Until now I was not able to find an obvious solution for this computations 
> using the FEEvaluation class. Maybe anyone else has a better idea for the 
> described problem.
>
> Thanks in advance,
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/7b3c1488-cf75-4526-8442-d1045dc73398n%40googlegroups.com.

[deal.II] Modifying shape function data in MatrixFree

['David' via deal.II User Group]

Hi all,

I'm currently trying to implement a vectorized variant of the residual 
assembly as it is done in step-44/one of the corresponding code-gallery 
examples using FEEvaluation objects in combination with matrix-free. 

I was not able to find an appropriate solution for the given code line 
and
 
the subsequent application for the assembly process: the major problem is 
that the shape function gradients are defined with regard to the current 
configuration, whereas my matrix-free object holds a mapping, which refers 
to the initial configuration. The reference configuration mapping  of my 
matrix-free object is desired as the final integration is actually 
performed in the reference frame.
A valid option is probably (I have not tried it) to use a matrix-free 
object with a different mapping (MappingQEulerian) which directly evaluates 
the shape gradients in the deformed configuration and use another 
referential matrix-free object for the integration part. 
The main reason to avoid this approach is that it would require a 
reinitialization of the 'deformed' matrix-free object in each nonlinear 
step and the mapping needs to be recomputed in each reinizialization. On 
the other hand, the expensive reinitialization is not really required as 
the mapping between the referential and the current configuration is known 
due to the (inverse) deformation gradient. 
Therefore, I'm looking for an opportunity to access the shape gradients and 
perform a push forward similarly to the approach in step-44 in order to 
evaluate the desired quantities. 

Until now I was not able to find an obvious solution for this computations 
using the FEEvaluation class. Maybe anyone else has a better idea for the 
described problem.

Thanks in advance,
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/34669b15-4c36-410c-9665-071b7bbe2c44n%40googlegroups.com.

[deal.II] Re: Solver for nonlinear solid dynamics

['David' via deal.II User Group]

Dear Wolfgang and Richard,

thanks a lot for your quick replies!

Have you taken a look at the deal.II Code Gallery? 

>https://dealii.org/developer/doxygen/deal.II/CodeGallery.html 
> There are a couple more large-deformation codes. 
>
>
Looked once, but obviously not close enough. There are several more 
appropriate codes I haven't noticed. 
 

> There is also this article from Thomas Wick about a deal.II-based code for 
> FSI: 
>https://media.archnumsoft.org/10305/ 
> I believe that the code that goes with the article is freely available as 
> well. 
>

Very nice! 



You might simplify one of the codes from the code gallery and plug in the 
> element integration routine from the last code to get to 'geometrically 
> nonlinear elasticity'.
> I have done some similar things, thus have a bit of code in that 
> direction. If you are interested in working together, just mail me.
>
> Good luck & Kind regards,
> Richard 
>

I think this is the way to go here. I would like to investigate this idea 
first of all and get familiar with the codes. Thanks for your kind offer, I 
will let you know, when it is needed.

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2d534a48-fc2e-4091-98b9-1b53271a88d8%40googlegroups.com.

[deal.II] Solver for nonlinear solid dynamics

['David' via deal.II User Group]

Hello everyone,

we (the preCICE team 
)
 
are targeting to support a new deal.II solver for preCICE. We want to 
simulate the probably well known Turek 

 
and Hron 
FSI3
 
benachmark 
.
 
Therefore, we would like to use a solver for nonlinear solid dynamics with 
deal.II , since the deformations are rather large.

This post is intended to collect some ideas before starting: Has anyone 
already developed such a solver or knows a good starting point?

We have currently a dynamic linear-elastic solid solver and there are of 
course several tutorial programs dealing with nonlinear static solid 
mechanics. 

Any idea is appreciated!

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/fd49415a-c9c1-4765-bdc5-45927ab55d9f%40googlegroups.com.