-
Notifications
You must be signed in to change notification settings - Fork 100
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Towards a generalized implementation of RT FEs that does not assume an oriented global mesh. #579
Towards a generalized implementation of RT FEs that does not assume an oriented global mesh. #579
Conversation
that does not assume an oriented global mesh. Not passing tests. To investigate.
Let us talk tomorrow at Monash. |
global DoFs are transformed into local DoFs
A follow-up on this PR. I have simplified the test above, using When I restrict the the interpolated function to cell number 2, I still get the function with the sign flipped, i.e., [-1 -2]. I have done all computations by hand (see attachment), and double-checked that they match what the program is doing. Thus, it might happen that I have missunderstood something in the implementation of RT FEs for arbitrarily oriented meshes @santiagobadia ... do you know what might be wrong in the notes? |
@santiagobadia ... Can it be that the signed determinant of the Jacobian has to be used for the Piola map if and only if you also use it when transforming the differential integral measure? I just thought that it makes sense as the Piola map somehow can be derived from the transformation of domain of integration. In any case, not sure about this. If this was true, it would mean that, e.g., fenics is not using the absolute value of the determinant of the Jacobian when transforming the Domain of integration .... How would you solve the missorientation issue when integrating? |
…g_rt_fes_for_arbitrarily_oriented_meshes
(documentation added to RaviartThomasRefFEs.jl) Two key ingredients in the implementation of this type of ReferenceFE are the get_shapefuns(model,cell_reffes) and get_dof_basis(mode,cell_reffes) overloads. These are written such that, for each cell K, they return the shape functions and dof values in the *global* RT space. For a dof owned by a face which is shared by two cells, there is a master and a slave cell. The slave cell first computes the shape functions and dof values using local-to-cell data structures, but then flips the sign of both in order to get their corresponding counterparts in the **global** RT space. As as result we have the following: * When we interpolate a function into the global FE space, and we perform the cell-wise DoF values to global DoF values gather operation, we can either extract the global DoF value from the master or slave cell without worrying about the sign. * When we evaluate a global FE function, and we perform the global DoF values to cell-wise DoF values scatter operation, we don't have to worry about the sign either. On the slave cell, we will have both the sign of the DoF value, and the sign of the shape function corresponding to the global DoF. * We do NOT have to use the signed determinant, but its absolute value, in the Piola Map.
@fverdugo @santiagobadia ... FYI ... I followed a different approach for generalizing RT FEs for arbitrarily oriented meshes to the one that I was trying to follow so far. See description of the following commit 3b6b880 (also added to the source code). I think that once the tests pass, the PR is ready to review. As a comment, I had to put the following code: in the |
Ok. I think this is no longer a problem. The code was using the signed determinant of the Jacobian both for the Piola Map and the calculation of the interior moments. Now we use the absolute value of the determinant both for the Piola Map and the calculation of interior moments. Thus, the outcome should be same, the tests with high order FEs should pass now. |
Codecov Report
@@ Coverage Diff @@
## master #579 +/- ##
==========================================
+ Coverage 87.55% 87.84% +0.29%
==========================================
Files 130 132 +2
Lines 13619 14183 +564
==========================================
+ Hits 11924 12459 +535
- Misses 1695 1724 +29
Continue to review full report at Codecov.
|
Ok, the tests passed. PR ready to review. |
1. We will assume that we will work with manifolds in which are cells are oriented as the global manifold. 2. For the Piola Map, cell_orientation is not actually necessary for RT FEs on Manifolds if want uses the approach in #579
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@amartinhuertas thanks for the PR!
I have some concerns about the introduced coupling between DiscreteModel and Raviart-Thomas reffes. See below.
src/FESpaces/ConformingFESpaces.jl
Outdated
ctype_reffe, cell_ctype = compress_cell_data(cell_reffe) | ||
ctype_num_dofs = map(num_dofs,ctype_reffe) | ||
ctype_ldof_comp = map(reffe->get_dof_to_comp(reffe),ctype_reffe) | ||
cell_conformity = CellConformity(cell_reffe) | ||
cell_shapefuns = lazy_map(get_shapefuns,cell_reffe,cell_map) | ||
cell_dof_basis = lazy_map(get_dof_basis,cell_reffe,cell_map) | ||
cell_shapefuns = get_shapefuns(model,cell_reffe) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you need to pass a full discrete model here? is not enough to do something like
cell_shapefuns = lazy_map(get_shapefuns,cell_reffe,cell_map, cell_flip)
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you need to pass a full discrete model here?
We shall agree that there is some work to do at the global FE space level to compute the (cell-wise restriction) of global shape functions and global DoFs in the case of global H(div)-conforming FE spaces (that, e.g., is not required to be done for global grad-conforming FE spaces). In the case of your code snippet, it would be to build the cell_flip
array out of the model. In the case of the code I have written, those lines of code which I added to ConformingFESpaces.jl
. The code which is in charge of doing this work cannot go to the ReferenceFEs
module, otherwise we would have a cycle. I guess that you are trying to avoid this cycle by having interfaces that do not depend on model, and that thus can go to ReferenceFEs.jl
. But I am not sure if this is a good solution, because of reasons which I try to make clear below.
I understand the concern of having code in ConformingFESpaces.jl
which depends on the particular kind of global FE space at hand. In order to solve this, I propose to add a new source file, DivConformingFEspaces.jl
, in FESpaces. See
https://github.com/gridap/Gridap.jl/pull/579/files#diff-ea5ed427dd3dd1fad920c1f2014e8950d5866705bde29f7278c6d93390ae72c1R1 for a motivation. In a nutshell, this source file implements those customizations of get_shape_funs
and get_dof_basis
for the global conforming FE space at hand.
is not enough to do something like
cell_shapefuns = lazy_map(get_shapefuns,cell_reffe,cell_map, cell_flip)
?
I am not sure if this would be convenient.
First, you would be passing the cell_flip
array in order to build the (cell-wise restriction) of global shape functions for all global conforming FE spaces, even if it is not needed (e.g., for Grad-Confirming spaces).
Second, it is hard to anticipate if cell_flip
is the info which will be needed in order to compute the (cell-wise restriction) of global shape functions for those global conforming FE spaces which we may want to implement in the future.
By having a get_shapefuns
interface with model, and the ReferenceFE, we have more flexibility in defining the global shape functions according to the global FE space at hand. If the name is confusing, we can perhaps think of another name, instead of reusing the one that is used in order to compute the reference FE shape functions.
Let me know what do you think.
src/FESpaces/ConformingFESpaces.jl
Outdated
get_cell_map(Triangulation(model))) | ||
end | ||
|
||
function get_dof_basis(model::DiscreteModel, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This includes a dependency between discrete model and RaviartThomas reffes. As a result, RaviartThomas reffes cannot be fully implemented in the ReferenceFEs module, which is sort of confusing. It is also a sign that the current API of ReferenceFEs is not general enough.
Do you think it's possible to fix this? Perhaps it is possible with some refactoring of the ReferenceFEs API.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See #579 (comment)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have been talking to @amartinhuertas this morning, and I think this is the way to go. Div-conforming FE spaces require to modify the reference basis when transformed to the local FE basis. That is the identity for Lagrangian FEM (grad-conforming) but not the case (in general) for other FEs like div or curl-conforming ones.
This is something that cannot be at the ref FE level (it is a global thing, to enforce div-conformity in the global space).
I also don't like to have dispatching wrt RT in the ConformingFESpace
file. We have also learned that the only modification that is required to go to div-conforming are these two functions that compute the global DOF and shape functions on physical cells. So, it seems to me that the most neat solution is to create a file for these FE spaces that simply implement these functions.
It is in the same spirit as e.g. discontinuous Galerkin and conforming FE spaces, which fill call the unconstrained FE constructor in different ways. We could also do the dispatching at that level for div-conforming, but since the code is the same with only the previous exception, I would just implement the functions that need dispatching.
Cheers,
--Santi
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We can deal with this If it is necessary and I would also rename the versions taking the discrete model to get_cell_shapefuns
and get_cell_dof_basis
like other functions returning cell-wise values.
Back in the day, get_shapefuns
and get_dof_basis
had only a single argument (the reference FE), but we included a second argument (the geometrical map) to handle RT and Nedelec. Now these FE types will be handled differently, so perhaps we can revert to the original definition of get_shapefuns
and get_dof_basis
with just a single argument. This would simplify the implementation of Lagrangian reffes, which became a bit dirty after introducing RT and Nedelec.
In summary, these would the API
# Abstract methods in ReferenceFEs module
get_shapefuns(reffe::ReferenceFE)
get_dof_basis(reffe::ReferenceFE)
# Methods with a default implementation in terms of previous ones (which will work for e.g. Lagrangian),
# but they can be overloaded by new ReferenceFE types
get_cell_shapefuns(model::DiscreteModel,cell_reffe::AbstractArray{<:ReferenceFE}))
get_cell_dof_basis(model::DiscreteModel,cell_reffe::AbstractArray{<:ReferenceFE}))
# Methods to be added for RT FEs
# 1) Helper methods specific for RT
_transformed_rt_shapefuns(reffe::GenericRefFE{RaviartThomas},map::Field,signflip::AbstractVector{<:Integer})
_transformed_rt_dof_basis(reffe::GenericRefFE{RaviartThomas},map::Field,signflip::AbstractVector{<:Integer})
# 2) Specialization of Abstract methods in ReferenceFEs module
get_shapefuns(reffe::GenericRefFE) # Already implemented
get_dof_basis(reffe::GenericRefFE) # Already implemented
# 3) Specialization of cell-wise (like already implemented in the PR)
get_cell_shapefuns(model::DiscreteModel,cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}}))
get_cell_dof_basis(model::DiscreteModel,cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}}))
What do you think?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Agreed
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Btw, the sign is not +1 or -1, not an Int
but an array with as many +1/-1 as DOFs in the RefFE
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok! I have updated the code snipped
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I used Bool
, not Int
. (true-> -1 false -> 1; not a big deal)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done in 1685e4d
Div-conforming FE spaces to a new source file
|
||
function get_sign_flip(model::DiscreteModel, | ||
cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}}) | ||
lazy_map(get_sign_flip, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
(Minor comment). Perhaps, it is easier to pre compute the the cell-wise sign flip into a standard julia vector of Int8
instead to try to make it lazy. It would be also significantly faster if you need to perform several loops over the vector, e.g. in non-linear problems
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The array is not as simple, it is an array of arrays of Int
. So, I would not proceed this way.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Then better lazy!
I would slightly refactor the definition of the kernel from which the lazy array is defined to match other similar kernels in Gridap. Here, the cell operation depends on a global state (the model). In this case, we usualy include the global state in a struct:
struct SignFlipMap{T} <: Map
model::T
end
function return_cache(k::SignFlipMap,reffe,cell)
model = k.model
# same code as below
end
function evaluate!(cache,k::SignFlipMap,reffe,cell)
# same code as below
end
# usage
lazy_map(SignFlipMap(model),cell_reffe,cell_ids)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done in 1685e4d
sign_flip = cache | ||
|
||
setsize!(sign_flip, (num_dofs(reffe),)) | ||
sign_flip .= false |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Replace this:
sign_flip = cache
setsize!(sign_flip, (num_dofs(reffe),))
by
setsize!(cache, (num_dofs(reffe),))
sign_flip = cache.array
Then sign_flip
(and the return of this function) will be a intrinsic julia array. This is also a common pattern used in Gridap to hide the cached vector from the user of this function.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done in 1685e4d
…g_rt_fes_for_arbitrarily_oriented_meshes
spaces \# Abstract methods in ReferenceFEs module get_shapefuns(reffe::ReferenceFE) get_dof_basis(reffe::ReferenceFE) \# Methods with a default implementation in terms of previous ones (which will work for e.g. Lagrangian), \# but they can be overloaded by new ReferenceFE types get_cell_shapefuns(model::DiscreteModel,cell_reffe::AbstractArray{<:ReferenceFE})) get_cell_dof_basis(model::DiscreteModel,cell_reffe::AbstractArray{<:ReferenceFE})) \# Methods to be added for global RT FE spaces \# 1) Helper methods specific for RT _transformed_rt_shapefuns(reffe::GenericRefFE{RaviartThomas},map::Field,signflip::AbstractVector{<:Integer}) _transformed_rt_dof_basis(reffe::GenericRefFE{RaviartThomas},map::Field,signflip::AbstractVector{<:Integer}) \# 2) Specialization of cell-wise (like already implemented in the PR) get_cell_shapefuns(model::DiscreteModel,cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}})) get_cell_dof_basis(model::DiscreteModel,cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}}))
@amartinhuertas the renaming It should be something like
There are the current functions in FESpace that are called |
After the renaming, you can introduce the new functions So the final API would be: get_fe_basis(fe_space) -> FEBasis
get_trial_fe_basis(fe_space) -> FEBasis
get_fe_dof_basis(fe_space) -> CellDof
get_cell_shapefuns(model,cell_reffes) -> AbstractArray{<:AbstractVector{<:Field}}
get_cell_dof_basis(model,cell_reffes) -> AbstractArray{<:AbstractVector{<:Dof}} |
…79/checks?check_run_id=2463353534#step:6:159 __collect_cell_space was calling the old interface of get_shapefuns
Ooops! You are right (too much rushing). Now it should be solved. Note also the following change that I had to do in order to pass the tests: I am not sure if it still makes sense to have |
Hi @amartinhuertas! The purpose of Why do you need to add |
Ok, thanks for the clarification. I think now I understand what you say. Thus, as I suspected (but didn't see), you were foreseeing uses cases that I am preventing with the last changes that I introduced (and should thus rollback).
Point 1. i in #579 (comment). See also #579 (comment) |
I would provide 2 constructors:
(in the 2 cases the second object can be optional since a default value can be taken from the first one) If the user wants to construct a FE space from a vector of |
One can even add these code into a 3rd constructor if necessary. |
Ok, I see. I will come back to you soon with (hopefully) the final version. |
a5f3d7f. Finally, we have two constructors of FESpace FESpace(model,cell_reffes,conformity=nothing) FESpace(model,cell_fe,cell_conformity::CellConformity) On the other hand, CellFE constructor becomes CellFE(model,cell_reffes,conformity=nothing)
Ok, done. Remarks:
|
The input should be a |
It is not actually needed. I can (and I will) definitely use a Conformity instance as a third parameter of the CellFE constructor. I did not do it because of two reasons:
In any case, I agree that using a Conformity parameter opens the door for (unforeseen) possibilities in the future, so that I will follow your advice. |
its third parameter
Done in 5e1afaf. Please continue your code review taking #579 (comment) as a basis. You will see there some questions that require an answer. |
src/FESpaces/FESpaceFactories.jl
Outdated
cell_fe::CellFE; | ||
conformity=nothing, | ||
cell_fe::CellFE, | ||
cell_conformity::CellConformity=CellConformity(cell_fe); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would not include a positional argument here. CellConformity
is in CellFE
and I don't think we want/consider the case in which we have different CellConformity
in CellFE
than the one passed in this constructor.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok, agreed. I do not see either why this is necessary, and can be even counterproductive.
Ok, done. I think the PR is more or less ready to be merged up to following comment that still needs to be addressed:
what do u think @fverdugo? |
I am in particular referring to these three methods: function CellConformity(cell_fe::CellFE,cell_conf::Nothing)
cell_fe.cell_conformity
end
function CellConformity(cell_fe::CellFE,cell_conf::CellConformity)
@assert length(cell_fe.cell_ctype) == length(cell_conf.cell_ctype)
cell_conf
end
# function CellConformity(cell_fe::CellFE,conformity::Union{Symbol,Conformity})
# if conformity in (L2Conformity(),:L2)
# @notimplemented
# elseif conformity == :default
# cell_fe.cell_conformity
# else
# @unreachable """\n
# Argument conformity should be either :default, :L2, or L2Conformity()
# """
# end
# end |
I agree that they are not longer needed. They can be removed. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi @amartinhuertas !
Some more comments!
I found that the name get_cell_shapefuns
was already correct in the Geometry
module (the new name get_fe_basis
does not make sense in this context).
See inlined-comments.
# shape function corresponding to the global DoF. | ||
# * We do NOT have to use the signed determinant, but its absolute value, in the Piola Map. | ||
|
||
function get_cell_dof_basis(model::DiscreteModel, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It looks that when using RT FEs with L2 conformity, the default cell shape funs will be generated (i.e. the ones without the transformation). Is this enough?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes this is enough because no conformity is to be imposed
@@ -70,9 +70,9 @@ end | |||
# lazy_append(a,b) | |||
#end | |||
|
|||
function get_cell_shapefuns(trian::AppendedTriangulation) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Here (in the Geometry
mdoule) the name get_cell_shapefuns
was already correct, since the result is the array of cell-wise shape functions.
Revert the renaming in this module and add import Gridap.Geometry: get_cell_shapefuns
in FESpaces.jl
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok, good call. Thanks for detecting this out!
test/CellDataTests/CellDofsTests.jl
Outdated
@@ -17,7 +17,7 @@ model = simplexify(CartesianDiscreteModel(domain,cells)) | |||
|
|||
trian = Triangulation(model) | |||
|
|||
v = GenericCellField(get_cell_shapefuns(trian),trian,ReferenceDomain()) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Idem as before, here the old name was correct!
get_fe_basis in the Geometry module
Done, ready to be accepted! (once the tests are passing, I will post a comment when that happens) |
All tests passing! Thanks @fverdugo and @santiagobadia for your support! I think we can now accept the PR. |
Pull request that stemmed from the discussion in #577
Tests are not passing. Spent quite a lot of time debugging the following simple case, and I do not see where is the source of the error.
@santiagobadia ... we should perhaps discuss what I did in this PR so that we can try to find the cause of the error.
In any case, I still do not understand two things: