div(v)*p and q*div(u) trick in reference space high level API for high order geometric mappings #636
Conversation
Codecov Report
@@ Coverage Diff @@
## master #636 +/- ##
==========================================
+ Coverage 87.66% 87.67% +0.01%
==========================================
Files 134 134
Lines 14267 14307 +40
==========================================
+ Hits 12507 12544 +37
- Misses 1760 1763 +3
Continue to review full report at Codecov.
|
|
Which is the math definition of the new operator DIV? It is more than the div wrt reference coordinates right? |
|
There is one paper by botchev ... Rehabilitation of lowest order ... That defines a discrete variant of the physical div as the normal flux integrated over the boundary divided by the element measure. I am not sure of the connections of that DIV operator and the one at hand (divergence in reference space) |
|
I want to understand why this new feature does not fit well in our current design. Perhaps it is related with the fact that the high level API always refers to the physical space (i.e. the gradient and div are always wrt the physical coordinate), but in some situations you also want to write integrals of gradients/divs in the reference cell with the high level API (which is not supported now) If this is the reason, perhaps we need to introduce a dΩ = Measure(Ω) # PhysicalDomain() by default
dω = Measure(Ω,ReferenceDomain())
∫( ∇(u)⋅∇(v) )dΩ # Gradients wrt physical coordinate
∫( ∇(u)⋅∇(v) )dω # Gradients wrt reference coordinateI like this API, but not sure how difficult it is to implement it. |
Ok! So the Which is the "tricky" part of the PR then? The |
|
Yes the reason you say. The high level API is designed to always perform the change of domain from physical to reference space. It computes the Jacobian of the mapping from reference to physical and uses it to transform the differential integral measure |
|
OK! then the part we have to work a bit more is the |
|
@fverdugo is not a trick, it is the right way to do it, there is a mathematically sound reason for all this. I like how it is implemented, but as you say, if you think there is a better way to deal with the What is the problem with the current design of this triangulation? |
|
I called it trick because by exploiting the particular definition of the pull back of RT shape function it lets you compute that integral without having to compute the divergence in physical space which is complicated for non affine geometries due to the piola map |
* Added a new trait template parameter to CellQuadrature (of type DomainStyle) to let one select whether to integrate in ReferenceDomain() or PhysicalDomain()
Ok, I already removed a6cc908#diff-37219e6468e24d7c4b29903834b894aa426a861aa49ef64122068cbd47ceb5daR20 BTW, another difference (benefit?) in the proposed solution is that we don't have different entries in the dictionary within |
There was a problem hiding this comment.
Hi @amartinhuertas thanks for the changes.
Just a very minor comment and it's ready to merge!
| cell_quad::AbstractArray{<:Quadrature} | ||
| cell_point::AbstractArray{<:AbstractArray{<:Point}} | ||
| cell_weight::AbstractArray{<:AbstractArray{<:Real}} | ||
| trian::Triangulation | ||
| domain_style::DS | ||
| data_domain_style::DDS | ||
| integration_domain_style::IDS | ||
| end |
There was a problem hiding this comment.
The default constructor of the struct has changed. Can you add an extra constructor with the old signature? and perhaps calling it in the tests? Just to be backwards compatible, I.e.,
function CellQuadrature(
cell_quad::AbstractArray{<:Quadrature},
cell_point::AbstractArray{<:AbstractArray{<:Point}},
cell_weight::AbstractArray{<:AbstractArray{<:Real}},
trian::Triangulation,
data_domain_style::DomainStyle)
CellQuadrature(cell_quad,cell_point,cell_weigth,trian,data_domain_style,PhysicalDomain())
endThere was a problem hiding this comment.
Yes! I will do that soon (I am in the middle of another task).
I guess that the current set of tests do not fail because they are not calling to the (old) default constructtor, right?
There was a problem hiding this comment.
Yes! I will do that soon (I am in the middle of another task).
Done!
Yes.... |
old CellQuadrature inner constructor for backward compatibility Added missing test
|
@fverdugo ... the PR can be accepted, right? |
|
Thanks! |
- Some fixes to DIV pending from PR #636
Hi @fverdugo,
I have come up with a proposal to support the following trick for the implementation of this term (Darcy+RT FEs)
This will let us by-pass the current limitation that we have in the computation of the divergence operator. (i..e, Gridap assumes that the mapping is affine, and thus that the Jacobian is constant within each cell). This prevents the usage of, e.g., quadrilateral bilinear disctrizations of the sphere using the so-called cubed-sphere mesh, but even worse, it prevents curved elements at all.
In code, you may se the new term in action here: b989735#diff-37219e6468e24d7c4b29903834b894aa426a861aa49ef64122068cbd47ceb5daR24
Concerns to take into account when reviewing:
DIV()operator is currently only defined for instances of typeSingleFieldFEBasis. Should it be defined for any other objects?ReferenceTriangulationtriangulation is implemented?ReferenceTriangulationTriangulationStylehere correct?TO-DO:
NEWS.mdonce we agree on the approach and the details.PS: @santiagobadia already took a look at the approach to implementation, he thinks it is reasonable, and agreed in the details.