Dear Martin,
I would like to start by thanking you for all of your help. Here is a link
to my fork: https://github.com/AlexanderCicchino/dealii
It is not fully working yet, I need your help for the "evaluate" as you
mentioned earlier. I have been trying multiple different ways but cannot
seem to get it to work properly.
First, I have a test setup in:
tests/mappings/mapping_q_generic_GCL_curvilinear.cc
where currently I am outputting everything, this will change in the future
and I will later on add a more complicated curvilinear mesh to test. But,
importantly, the current implementation fails the test.
Also, I made changes to source/fe/mapping_q_generic.cc as you suggested in
the function mentioned above. The last part, which is including the
evaluate is in:
https://github.com/AlexanderCicchino/dealii/blob/master/source/fe/mapping_q_generic.cc
line 1622. We have X_l*grad(X_m) evaluated at each quadrature point, now we
need to somehow have the gradient of that evaluated at the quadrature
points. I assumed after line 1622 that that gradient is written in
"grad_Xl_grad_Xm" which is gradient( X_l * gradient(X_m) ) then I loop
cyclically through it for the conservative curl form.
Please let me know how you suggest I should proceed/setup the evaluate call
at line 1624.
Also, I noticed that the conservative curl form from Kopriva is not well
posed for 2D. In the past, for 2D we would extend the grid by unit 1 in the
z direction to properly evaluate the metric terms since, for example we
need: d/(d \zeta) ( z* dy / (d\eta) ). Any suggestions on how to implement
the 2D version in mapping Q generic?
Thank you,
Alex
On Wednesday, June 24, 2020 at 9:40:55 AM UTC-4, Martin Kronbichler wrote:
>
> Dear Alex,
>
> Great! I would suggest to start by simply adding new code to the
> maybe_update_q_points_Jacobians_... function with the option to turn it off
> or on. Depending on how the final implementation will look like we might
> want to move that to a separate place, but I think it will be less
> repetitive if we use a single place.
>
> Best,
> Martin
> On 22.06.20 19:59, Alexander Cicchino wrote:
>
> Dear Martin,
>
> Thank you very much! I have been working on making the test case not
> depend on our in house flowsolver's functions.
> I think that implementing Eq. 36 the "conservative curl" form would be
> sufficient.
> Yes this procedure sounds perfect to me, and I agree with the dimension of
> the object described. I have been going through the source code that you
> sent to familiarize myself with the objects. Should I be adding to the
> function maybe_update_q_points_Jacobians_and_grads_tensor or should I
> create a new function for it?
>
> Thank you,
> Alex
>
> On Friday, June 19, 2020 at 5:09:14 AM UTC-4, Martin Kronbichler wrote:
>>
>> Dear Alex,
>>
>> Great! The first thing we need to know is the equation. I had a quick
>> look in the paper by Kopriva and I think we want to use either equation
>> (36) or (37), depending on whether we consider the conservative or
>> invariant curl form, respectively. In either case, it appears that we need
>> to do this in a two-step procedure. The first step is to compute X_l and
>> \nabla_\xi X_m, which in deal.II speak are the "q_points" and "Jacobians".
>> The implementation in mapping_q_generic.cc is a bit involved because we
>> have a slow algorithm (working for arbitrary quadrature rules) and a fast
>> one for tensor product quadrature rules. Let us consider the fast one
>> because I think we have most ingredients available, whereas we would need
>> to fill additional fields for the slow algorithm. The source code for those
>> parts is here:
>>
>>
>> https://github.com/dealii/dealii/blob/9e05a87db802ecd073bf7567d77f3491170d84b4/source/fe/mapping_q_generic.cc#L1463-L1592
>>
>> I skipped the part on the Hessians (second derivative of transformation)
>> because we won't need them. The important parts here are the extractions of
>> the positions in line 1554 and the extraction of the gradients
>> (contravariant transformations) in line 1575. Those two parts will be the
>> starting point for the second phase we need to do in addition: According to
>> the algorithm by Kopriva, we need to define this as the curl of the
>> discrete interpolation of X_l \nabla_\xi X_m. To get the curl, we need
>> another round through the SelectEvaluator::evaluate() call in that function
>> to get the reference-cell gradient of that object, from which we can then
>> collect the entries of the curl. To call into evaluate one more time, we
>> also need a new data.shape_info object that does the collocation evaluation
>> of derivatives. That should only be two lines that I can show you how and
>> where to add, so let us not worry about that part now. What is important to
>> understand first (in terms of index notation vs tensor notation) is the
>> dimension of the object. I believe that X_l \nabla_\xi X_m is a rank-two
>> tensor, so it has
