Re: [deal.II] Geometric Conservation Law

[alexande...@gmail.com]

Dear all,
To complete this topic for both MappingQGeneric and MappingFEField, I am 
wondering if it is possible to extract the quadrature in the reference 
element in the function maybe_update_Jacobians and maybe_compute_q_points? 
These are needed to ensure that the normals are consistent at each surface 
quadrature node with the metric terms that satisfy GCL. 
I was able to get the reference quadrature nodes in 
maybe_update_q_points_Jacobians_and_grads_tensor in fe/mapping_q_generic.cc 
by using data.shape_info.data[0].quadrature, but when it is not a tensor 
product basis (in the 2 functions listed above for mapping_q_generic.cc) I 
noticed that data.shape_info is empty (and there is no data.shape_info in 
mapping_fe_field.cc). Ideally, is it possible to fill data.shape_info in 
mapping_q_generic.cc or are there any other suggestions on how to extract 
the quadrature on the reference element for both mapping_q_generic.cc and 
mapping_fe_field.cc (just the 1D representation is sufficient)?
Thank you,
Alex


On Friday, July 10, 2020 at 12:37:21 PM UTC-4 alexande...@gmail.com wrote:

> Dear Martin or whom it may concern,
>
> I have solved the previous problem and am confirming that the conservative 
> curl form has now been implemented and passes 2 complicated tests for GCL 
> on symmetric and non-symmetric curvilinear grids for different polynomial 
> degrees. Turns out that we did not need to use the "evaluate" call. The 
> changes are found in:
>
> https://github.com/AlexanderCicchino/dealii/blob/master/source/fe/mapping_q_generic.cc
> lines 1599 - 1697 (line 1602 to switch it on or off)
> and the 2 tests are /tests/mappings/mapping_q_generic_GCL_curvilinear.cc 
> for symmetric curv. grid and 
> /tests/mappings/mapping_q_generic_GCL_curvilinear_nonsym.cc for non-sym. 
> curv grid.
> How should I proceed with the pull request?
>
> Also, for additional work, are there any suggestions on how to proceed 
> with 2D. Additionally, are there any suggestions for templating the mapping 
> for ease automatically differentiating the metric terms?
>
> Thank you for all the help,
> Alex
>
>
> On Thursday, July 9, 2020 at 11:54:12 AM UTC-4, Alexander Cicchino wrote:
>>
>> Dear Martin,
>>
>> As an update, I figured out a way to have it work in 3D without needing 
>> to use another Evaluate call. The only issue is that I need the quadrature 
>> in the reference element since I construct an FE_DQGArbitraryNodes to 
>> evaluate the derivative (I do this since the metric Jacobian is in a sense 
>> "collocated" on the quadrature points so the gradient needed is just the 
>> derivative of the basis functions). 
>> I have it passing collocated and uncollocated curvilinear test cases in 
>> 3D now perfectly, but I currently construct the FE_DGQArbitraryNodes in 
>> mapping_q_generic.cc by explicitly constructing a Quadrature like the 
>> testcase (example QGauss if the testcase uses QGauss). 
>>
>> I couldn't find a way to extract the reference element quadrature from 
>> the MappingQGeneric InternalData, do you have any suggestions?
>>
>> Thank you,
>> Alex
>>
>> On Friday, July 3, 2020 at 6:11:28 PM UTC-4, Alexander Cicchino wrote:
>>>
>>> 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,

Re: [deal.II] Geometric Conservation Law

[Alexander Cicchino]

Dear Martin or whom it may concern,

I have solved the previous problem and am confirming that the conservative 
curl form has now been implemented and passes 2 complicated tests for GCL 
on symmetric and non-symmetric curvilinear grids for different polynomial 
degrees. Turns out that we did not need to use the "evaluate" call. The 
changes are found in:
https://github.com/AlexanderCicchino/dealii/blob/master/source/fe/mapping_q_generic.cc
lines 1599 - 1697 (line 1602 to switch it on or off)
and the 2 tests are /tests/mappings/mapping_q_generic_GCL_curvilinear.cc 
for symmetric curv. grid and 
/tests/mappings/mapping_q_generic_GCL_curvilinear_nonsym.cc for non-sym. 
curv grid.
How should I proceed with the pull request?

Also, for additional work, are there any suggestions on how to proceed with 
2D. Additionally, are there any suggestions for templating the mapping for 
ease automatically differentiating the metric terms?

Thank you for all the help,
Alex


On Thursday, July 9, 2020 at 11:54:12 AM UTC-4, Alexander Cicchino wrote:
>
> Dear Martin,
>
> As an update, I figured out a way to have it work in 3D without needing to 
> use another Evaluate call. The only issue is that I need the quadrature in 
> the reference element since I construct an FE_DQGArbitraryNodes to evaluate 
> the derivative (I do this since the metric Jacobian is in a sense 
> "collocated" on the quadrature points so the gradient needed is just the 
> derivative of the basis functions). 
> I have it passing collocated and uncollocated curvilinear test cases in 3D 
> now perfectly, but I currently construct the FE_DGQArbitraryNodes in 
> mapping_q_generic.cc by explicitly constructing a Quadrature like the 
> testcase (example QGauss if the testcase uses QGauss). 
>
> I couldn't find a way to extract the reference element quadrature from the 
> MappingQGeneric InternalData, do you have any suggestions?
>
> Thank you,
> Alex
>
> On Friday, July 3, 2020 at 6:11:28 PM UTC-4, Alexander Cicchino wrote:
>>
>> 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 Ale

Re: [deal.II] Geometric Conservation Law

[Alexander Cicchino]

Dear Martin,

As an update, I figured out a way to have it work in 3D without needing to 
use another Evaluate call. The only issue is that I need the quadrature in 
the reference element since I construct an FE_DQGArbitraryNodes to evaluate 
the derivative (I do this since the metric Jacobian is in a sense 
"collocated" on the quadrature points so the gradient needed is just the 
derivative of the basis functions). 
I have it passing collocated and uncollocated curvilinear test cases in 3D 
now perfectly, but I currently construct the FE_DGQArbitraryNodes in 
mapping_q_generic.cc by explicitly constructing a Quadrature like the 
testcase (example QGauss if the testcase uses QGauss). 

I couldn't find a way to extract the reference element quadrature from the 
MappingQGeneric InternalData, do you have any suggestions?

Thank you,
Alex

On Friday, July 3, 2020 at 6:11:28 PM UTC-4, Alexander Cicchino wrote:
>
> 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

Re: [deal.II] Geometric Conservation Law

[Alexander Cicchino]

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

Re: [deal.II] Geometric Conservation Law

[Martin Kronbichler]

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 dim*dim components, and we compute the
gradient that gives us a dim * dim * dim tensor. Taking the curl
in the derivative and inner tensor dimension space, we get rid of
one component and up with a dim * dim tensor as expected. The
final step we need to do is to divide by the determinant of the
Jacobian (contravariant vectors), because the inverse Jacobian in
deal.II does not yet pre-multiply with the determinant.

Does that procedure sound reasonable to you? If yes, we could
start putting together the ingredients. It would be good to have a
mesh (the curvilinear case you were mentioning) where we can test
those formulas.

Best,
Martin

On 17.06.20 18:37, Alexander Cicchino wrote:

Dear Martin,

Thank you for your response. Yes I agree that only some local
computations are necessary to implement the identities.
Yes I would be interested in this feature and trying to implement
it. Do you have any suggestions on where I should start and
overall practices I should follow?

Thank you,
Alex

On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin
Kronbichler wrote:

Dear Alex,

This has been on my list of things to implement and verify
with deal.II over a range of examples for quite a while, so
I'm glad you bringing the topic up. It is definitely true
that our way to define Jacobians does not take those
identities into account, but I believe we should add support
for them. The nice thing is that only some local computations
are necessary, so having the option to use it in

Re: [deal.II] Geometric Conservation Law

[Alexander Cicchino]

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 dim*dim components, and we compute the gradient that 
> gives us a dim * dim * dim tensor. Taking the curl in the derivative and 
> inner tensor dimension space, we get rid of one component and up with a dim 
> * dim tensor as expected. The final step we need to do is to divide by the 
> determinant of the Jacobian (contravariant vectors), because the inverse 
> Jacobian in deal.II does not yet pre-multiply with the determinant.
>
> Does that procedure sound reasonable to you? If yes, we could start 
> putting together the ingredients. It would be good to have a mesh (the 
> curvilinear case you were mentioning) where we can test those formulas.
>
> Best,
> Martin
> On 17.06.20 18:37, Alexander Cicchino wrote:
>
> Dear Martin,
>
> Thank you for your response. Yes I agree that only some local computations 
> are necessary to implement the identities.
> Yes I would be interested in this feature and trying to implement it. Do 
> you have any suggestions on where I should start and overall practices I 
> should follow?
>
> Thank you,
> Alex
>
> On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin Kronbichler wrote: 
>>
>> Dear Alex,
>>
>> This has been on my list of things to implement and verify with deal.II 
>> over a range of examples for quite a while, so I'm glad you bringing the 
>> topic up. It is definitely true that our way to define Jacobians does not 
>> take those identities into account, but I believe we should add support for 
>> them. The nice thing is that only some local computations are necessary, so 
>> having the option to use it in the polynomial mapping classes would be 
>> great. If you would be interested in this feature and trying to implement 
>> things, I'd be happy to guide you to the right places in the code.
>>
>> Best,
>> Martin
>> On 17.06.20 06:02, Alexander Cicchino wrote:
>>
>> Thank you for responding Wolfgang Bangerth.
>>
>> The GCL condition comes from the discretized scheme satisfying 
>> free-stream preservation. I will demonstrate this for 2D below, (can be 
>> interpreted for spectral, DG, finite difference, finite volume etc):
>> Consider the conservation law: \frac{\partial W}{\partial t} + 
>> \frac{\partial F}{\partial x} +\frac{\p

Re: [deal.II] Geometric Conservation Law

[Martin Kronbichler]

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 dim*dim 
components, and we compute the gradient that gives us a dim * dim * dim 
tensor. Taking the curl in the derivative and inner tensor dimension 
space, we get rid of one component and up with a dim * dim tensor as 
expected. The final step we need to do is to divide by the determinant 
of the Jacobian (contravariant vectors), because the inverse Jacobian in 
deal.II does not yet pre-multiply with the determinant.


Does that procedure sound reasonable to you? If yes, we could start 
putting together the ingredients. It would be good to have a mesh (the 
curvilinear case you were mentioning) where we can test those formulas.


Best,
Martin

On 17.06.20 18:37, Alexander Cicchino wrote:

Dear Martin,

Thank you for your response. Yes I agree that only some local 
computations are necessary to implement the identities.
Yes I would be interested in this feature and trying to implement it. 
Do you have any suggestions on where I should start and overall 
practices I should follow?


Thank you,
Alex

On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin Kronbichler 
wrote:


Dear Alex,

This has been on my list of things to implement and verify with
deal.II over a range of examples for quite a while, so I'm glad
you bringing the topic up. It is definitely true that our way to
define Jacobians does not take those identities into account, but
I believe we should add support for them. The nice thing is that
only some local computations are necessary, so having the option
to use it in the polynomial mapping classes would be great. If you
would be interested in this feature and trying to implement
things, I'd be happy to guide you to the right places in the code.

Best,
Martin

On 17.06.20 06:02, Alexander Cicchino wrote:

Thank you for responding Wolfgang Bangerth.

The GCL condition comes from the discretized scheme satisfying
free-stream preservation. I will demonstrate this for 2D below,
(can be interpreted for spectral, DG, finite difference, finite
volume etc):
Consider the conservation law: \frac{\partial W}{\partial t} +
\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0
Transforming this to the reference computational space
(x,y)->(\xi, \eta):
J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial
x} * \frac{\partial F}{\partial \xi} + J * \frac{ \partial
\eta}{\partial x}* \frac{\partial F}{\partial \eta} + J * \frac{
\partial \xi}{\partial y} * \frac{\partial G}{\partial \xi} +
J*\frac{ \partial \eta}{\partial y}*\frac{\partial G}{\partial \eta}
Putting this in conservative form results in:
J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} (
J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial
\xi}{\partial y} ) + \frac{\partial}{\partial \eta} (
J*F*\frac{\partial \eta}{\partial x} +

Re: [deal.II] Geometric Conservation Law

[Alexander Cicchino]

Dear Martin,

Thank you for your response. Yes I agree that only some local computations 
are necessary to implement the identities.
Yes I would be interested in this feature and trying to implement it. Do 
you have any suggestions on where I should start and overall practices I 
should follow?

Thank you,
Alex

On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin Kronbichler wrote:
>
> Dear Alex,
>
> This has been on my list of things to implement and verify with deal.II 
> over a range of examples for quite a while, so I'm glad you bringing the 
> topic up. It is definitely true that our way to define Jacobians does not 
> take those identities into account, but I believe we should add support for 
> them. The nice thing is that only some local computations are necessary, so 
> having the option to use it in the polynomial mapping classes would be 
> great. If you would be interested in this feature and trying to implement 
> things, I'd be happy to guide you to the right places in the code.
>
> Best,
> Martin
> On 17.06.20 06:02, Alexander Cicchino wrote:
>
> Thank you for responding Wolfgang Bangerth.
>
> The GCL condition comes from the discretized scheme satisfying free-stream 
> preservation. I will demonstrate this for 2D below, (can be interpreted for 
> spectral, DG, finite difference, finite volume etc):
> Consider the conservation law: \frac{\partial W}{\partial t} + 
> \frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0
> Transforming this to the reference computational space (x,y)->(\xi, \eta):
> J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial x} * 
> \frac{\partial F}{\partial \xi} + J * \frac{ \partial \eta}{\partial x}* 
> \frac{\partial F}{\partial \eta} + J * \frac{ \partial \xi}{\partial y} * 
> \frac{\partial G}{\partial \xi} + J*\frac{ \partial \eta}{\partial 
> y}*\frac{\partial G}{\partial \eta}
> Putting this in conservative form results in:
> J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} ( 
> J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial \xi}{\partial y} ) 
> + \frac{\partial}{\partial \eta} ( J*F*\frac{\partial \eta}{\partial x} 
> +J*G*\frac{\partial \eta}{\partial y} ) - F*( GCL in x) - G*(GCL in y) =0
>
> where GCL in x = \frac{\partial }{\partial \xi} ( det(J)* \frac{\partial 
> \xi 
>  }{\partial x}) + \frac{\partial }{\partial \eta}( det(J)* \frac{\partial 
> \eta}{\partial x} )
> similarly for y.
>
> So for the conservative numerical scheme to satisfy free stream 
> preservation, the GCL conditions must go to zero.
> For linear grids, there are no issues with the classical definition for 
> the inverse of the Jacobian, but what Kopriva had shown (before him Thomas 
> and Lombard), was that the metric Jacobian has to be calculated in either a 
> "conservative curl form" or an "invariant curl form" since it reduces the 
> GCL condition to the divergence of a curl, which is always discretely 
> satisfied. In the paper by Kopriva, he shows this, an example in 3D:
>  Analytically
> J*\frac{\partial \xi}{\partial x} = \frac{\partial z}{\partial \zeta} * 
> \frac{\partial y}{\partial \eta} - \frac{\partial z}{\partial \eta} * 
> \frac{\partial y}{\partial \zeta}
>
> but the primer doesn't satisfy free-stream preservation while the latter 
> ("conservative curl form") does.
>
> I will put together a unit test for a curvilinear grid. 
>
> Thank you,
> Alex
>
> On Tuesday, June 16, 2020 at 10:24:59 PM UTC-4, Wolfgang Bangerth wrote: 
>>
>>
>> Alexander, 
>>
>> > I am wondering if anybody has also found that the inverse of the 
>> Jacobian from 
>> > FE Values, with MappingQGeneric does not satisfy the Geometric 
>> Conservation 
>> > Law (GCL), in the sense of: 
>> > 
>> > Kopriva, David A. "Metric identities and the discontinuous spectral 
>> element 
>> > method on curvilinear meshes." /Journal of Scientific Computing/ 26.3 
>> (2006): 301. 
>> > 
>> > on curvilinear elements/manifolds in 3D. 
>> > That is: 
>> > \frac{\partial }{\partial \hat{x}_1} *det(J)* \frac{\partial \hat{x}_1 
>> > }{\partial x_1} + \frac{\partial }{\partial \hat{x}_2} *det(J)* 
>> \frac{\partial 
>> > \hat{x}_2}{\partial x} + \frac{\partial }{\partial \hat{x}_3} * 
>> > det(J)*\frac{\partial \hat{x}_3 }{\partial x_1} != 0 (GCL says it 
>> should =0, 
>> > similarly for x_2 and x_3) 
>> > 
>> > If so or if not, also, has anybody found a remedy to have the inverse 
>> of the 
>> > Jacobian from FE Values with MappingQGeneric to satisfy the GCL. 
>>
>> I'm not sure any of us have ever thought about it. (I haven't -- but I 
>> really 
>> shouldn't speak for anyone else.) Can you explain what this equality 
>> represents? Why should it hold? 
>>
>> I'm also unsure whether we've ever checked whether it holds (exactly or 
>> approximately). Can you create a small test program that illustrates the 
>> behavior you are seeing? 
>>
>> Best 
>>   W. 
>>
>> -- 
>>  
>> Wolfgang Bangerth  

Re: [deal.II] Geometric Conservation Law

[Martin Kronbichler]

Dear Alex,

This has been on my list of things to implement and verify with deal.II 
over a range of examples for quite a while, so I'm glad you bringing the 
topic up. It is definitely true that our way to define Jacobians does 
not take those identities into account, but I believe we should add 
support for them. The nice thing is that only some local computations 
are necessary, so having the option to use it in the polynomial mapping 
classes would be great. If you would be interested in this feature and 
trying to implement things, I'd be happy to guide you to the right 
places in the code.


Best,
Martin

On 17.06.20 06:02, Alexander Cicchino wrote:

Thank you for responding Wolfgang Bangerth.

The GCL condition comes from the discretized scheme satisfying 
free-stream preservation. I will demonstrate this for 2D below, (can 
be interpreted for spectral, DG, finite difference, finite volume etc):
Consider the conservation law: \frac{\partial W}{\partial t} + 
\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0

Transforming this to the reference computational space (x,y)->(\xi, \eta):
J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial x} * 
\frac{\partial F}{\partial \xi} + J * \frac{ \partial \eta}{\partial 
x}* \frac{\partial F}{\partial \eta} + J * \frac{ \partial 
\xi}{\partial y} * \frac{\partial G}{\partial \xi} + J*\frac{ \partial 
\eta}{\partial y}*\frac{\partial G}{\partial \eta}

Putting this in conservative form results in:
J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} ( 
J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial \xi}{\partial 
y} ) + \frac{\partial}{\partial \eta} ( J*F*\frac{\partial 
\eta}{\partial x} +J*G*\frac{\partial \eta}{\partial y} ) - F*( GCL in 
x) - G*(GCL in y) =0


where GCL in x = \frac{\partial }{\partial \xi} ( det(J)* 
\frac{\partial \xi

 }{\partial x}) + \frac{\partial }{\partial \eta}( det(J)* \frac{\partial
\eta}{\partial x} )
similarly for y.

So for the conservative numerical scheme to satisfy free stream 
preservation, the GCL conditions must go to zero.
For linear grids, there are no issues with the classical definition 
for the inverse of the Jacobian, but what Kopriva had shown (before 
him Thomas and Lombard), was that the metric Jacobian has to be 
calculated in either a "conservative curl form" or an "invariant curl 
form" since it reduces the GCL condition to the divergence of a curl, 
which is always discretely satisfied. In the paper by Kopriva, he 
shows this, an example in 3D:

 Analytically
J*\frac{\partial \xi}{\partial x} = \frac{\partial z}{\partial \zeta} 
* \frac{\partial y}{\partial \eta} - \frac{\partial z}{\partial \eta} 
* \frac{\partial y}{\partial \zeta}


but the primer doesn't satisfy free-stream preservation while the 
latter ("conservative curl form") does.


I will put together a unit test for a curvilinear grid.

Thank you,
Alex

On Tuesday, June 16, 2020 at 10:24:59 PM UTC-4, Wolfgang Bangerth wrote:


Alexander,

> I am wondering if anybody has also found that the inverse of the
Jacobian from
> FE Values, with MappingQGeneric does not satisfy the Geometric
Conservation
> Law (GCL), in the sense of:
>
> Kopriva, David A. "Metric identities and the discontinuous
spectral element
> method on curvilinear meshes." /Journal of Scientific Computing/
26.3 (2006): 301.
>
> on curvilinear elements/manifolds in 3D.
> That is:
> \frac{\partial }{\partial \hat{x}_1} *det(J)* \frac{\partial
\hat{x}_1
> }{\partial x_1} + \frac{\partial }{\partial \hat{x}_2} *det(J)*
\frac{\partial
> \hat{x}_2}{\partial x} + \frac{\partial }{\partial \hat{x}_3} *
> det(J)*\frac{\partial \hat{x}_3 }{\partial x_1} != 0 (GCL says
it should =0,
> similarly for x_2 and x_3)
>
> If so or if not, also, has anybody found a remedy to have the
inverse of the
> Jacobian from FE Values with MappingQGeneric to satisfy the GCL.

I'm not sure any of us have ever thought about it. (I haven't --
but I really
shouldn't speak for anyone else.) Can you explain what this equality
represents? Why should it hold?

I'm also unsure whether we've ever checked whether it holds
(exactly or
approximately). Can you create a small test program that
illustrates the
behavior you are seeing?

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

Re: [deal.II] Geometric Conservation Law

[Alexander Cicchino]

Thank you for responding Wolfgang Bangerth.

The GCL condition comes from the discretized scheme satisfying free-stream 
preservation. I will demonstrate this for 2D below, (can be interpreted for 
spectral, DG, finite difference, finite volume etc):
Consider the conservation law: \frac{\partial W}{\partial t} + 
\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0
Transforming this to the reference computational space (x,y)->(\xi, \eta):
J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial x} * 
\frac{\partial F}{\partial \xi} + J * \frac{ \partial \eta}{\partial x}* 
\frac{\partial F}{\partial \eta} + J * \frac{ \partial \xi}{\partial y} * 
\frac{\partial G}{\partial \xi} + J*\frac{ \partial \eta}{\partial 
y}*\frac{\partial G}{\partial \eta}
Putting this in conservative form results in:
J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} ( 
J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial \xi}{\partial y} ) 
+ \frac{\partial}{\partial \eta} ( J*F*\frac{\partial \eta}{\partial x} 
+J*G*\frac{\partial \eta}{\partial y} ) - F*( GCL in x) - G*(GCL in y) =0

where GCL in x = \frac{\partial }{\partial \xi} ( det(J)* \frac{\partial 
\xi 
 }{\partial x}) + \frac{\partial }{\partial \eta}( det(J)* \frac{\partial 
\eta}{\partial x} )
similarly for y.

So for the conservative numerical scheme to satisfy free stream 
preservation, the GCL conditions must go to zero.
For linear grids, there are no issues with the classical definition for the 
inverse of the Jacobian, but what Kopriva had shown (before him Thomas and 
Lombard), was that the metric Jacobian has to be calculated in either a 
"conservative curl form" or an "invariant curl form" since it reduces the 
GCL condition to the divergence of a curl, which is always discretely 
satisfied. In the paper by Kopriva, he shows this, an example in 3D:
 Analytically
J*\frac{\partial \xi}{\partial x} = \frac{\partial z}{\partial \zeta} * 
\frac{\partial y}{\partial \eta} - \frac{\partial z}{\partial \eta} * 
\frac{\partial y}{\partial \zeta}

but the primer doesn't satisfy free-stream preservation while the latter 
("conservative curl form") does.

I will put together a unit test for a curvilinear grid. 

Thank you,
Alex

On Tuesday, June 16, 2020 at 10:24:59 PM UTC-4, Wolfgang Bangerth wrote:
>
>
> Alexander, 
>
> > I am wondering if anybody has also found that the inverse of the 
> Jacobian from 
> > FE Values, with MappingQGeneric does not satisfy the Geometric 
> Conservation 
> > Law (GCL), in the sense of: 
> > 
> > Kopriva, David A. "Metric identities and the discontinuous spectral 
> element 
> > method on curvilinear meshes." /Journal of Scientific Computing/ 26.3 
> (2006): 301. 
> > 
> > on curvilinear elements/manifolds in 3D. 
> > That is: 
> > \frac{\partial }{\partial \hat{x}_1} *det(J)* \frac{\partial \hat{x}_1 
> > }{\partial x_1} + \frac{\partial }{\partial \hat{x}_2} *det(J)* 
> \frac{\partial 
> > \hat{x}_2}{\partial x} + \frac{\partial }{\partial \hat{x}_3} * 
> > det(J)*\frac{\partial \hat{x}_3 }{\partial x_1} != 0 (GCL says it should 
> =0, 
> > similarly for x_2 and x_3) 
> > 
> > If so or if not, also, has anybody found a remedy to have the inverse of 
> the 
> > Jacobian from FE Values with MappingQGeneric to satisfy the GCL. 
>
> I'm not sure any of us have ever thought about it. (I haven't -- but I 
> really 
> shouldn't speak for anyone else.) Can you explain what this equality 
> represents? Why should it hold? 
>
> I'm also unsure whether we've ever checked whether it holds (exactly or 
> approximately). Can you create a small test program that illustrates the 
> behavior you are seeing? 
>
> 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/4f313231-dbb3-445f-923c-9eaff17ab783o%40googlegroups.com.

Re: [deal.II] Geometric Conservation Law

[Wolfgang Bangerth]

Alexander,

I am wondering if anybody has also found that the inverse of the Jacobian from 
FE Values, with MappingQGeneric does not satisfy the Geometric Conservation 
Law (GCL), in the sense of:


Kopriva, David A. "Metric identities and the discontinuous spectral element 
method on curvilinear meshes." /Journal of Scientific Computing/ 26.3 (2006): 301.


on curvilinear elements/manifolds in 3D.
That is:
\frac{\partial }{\partial \hat{x}_1} *det(J)* \frac{\partial \hat{x}_1 
}{\partial x_1} + \frac{\partial }{\partial \hat{x}_2} *det(J)* \frac{\partial 
\hat{x}_2}{\partial x} + \frac{\partial }{\partial \hat{x}_3} * 
det(J)*\frac{\partial \hat{x}_3 }{\partial x_1} != 0 (GCL says it should =0, 
similarly for x_2 and x_3)


If so or if not, also, has anybody found a remedy to have the inverse of the 
Jacobian from FE Values with MappingQGeneric to satisfy the GCL.


I'm not sure any of us have ever thought about it. (I haven't -- but I really 
shouldn't speak for anyone else.) Can you explain what this equality 
represents? Why should it hold?


I'm also unsure whether we've ever checked whether it holds (exactly or 
approximately). Can you create a small test program that illustrates the 
behavior you are seeing?


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/d8421b9f-0fa9-9769-b63c-281500879687%40colostate.edu.