Skip to content
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

Misc changes in order to support RT FEs in general manifolds (i.e., num_cell_dims(model) < num_point_dims(model)) #577

Merged
merged 11 commits into from
May 7, 2021

Conversation

amartinhuertas
Copy link
Member

@amartinhuertas amartinhuertas commented Apr 11, 2021

Hi @santiagobadia, @fverdugo,

In this PR I include a first proposal for a set of miscellaneous changes/fixes required in order to support RT FEs posed over general manifolds that I want to be subject to review/discussion.

These changes are related, among others, to the remarks given in issue #571 and references therein.

As far as I understand, the model should eventually provide the cell_orientation array. We should also discuss how can we include this info in the abstract interface of DiscreteModel and related data types, and how this information is computed in the case of general, unstructured, oriented manifolds.

Apart from using cell_orientation for the implementation of the sign of the ContraVariant Piola Map, I am wondering where (if at any place) this info should also be used. For example, should it be used for the evaluation of dofs, and in particular, for the computation of normals in physical space from those of reference space? I dont think so. However, I am concerned that, in FENICS, this info seems to be required for the evaluation of dofs (see, e.g., https://fenics-dolfin.readthedocs.io/en/2017.2.0/apis/api_ufc.html#_CPPv2NK3ufc14finite_element3ufc14finite_element13evaluate_dofsEPdRK8functionPKdiRK4cellPKN3ufc18coordinate_mappingE) apart from a myriad of other subroutines. Am I missing something?

Currently I am assuming that all cells are oriented as the global manifold in the case of num_cell_dims(model) < num_point_dims(model), i.e., cell_orientation=false for all cells. In the case of num_cell_dims(model) == num_point_dims(model) I am extracting this info from the sign of the determinant of the Jacobian. Regarding this last point, although it is temporary, I am concerned about the fact that I believe that we are computing the Jacobian and the determinant for both the ContraVariantPiolaMap and cell_orientation. Is that true? If yes, is there anyway to waive this duplicated computation? On the other hand, the Jacobian and its (pseudo)determinant are also required for the transformation of the integral differential measure (apart from the Piola transformation). Does Gridap ensure that this is only computed once? If yes, which is the mechanism by which it is achieving this?

(i.e., num_cell_dims(model) < num_point_dims(model))
@amartinhuertas
Copy link
Member Author

Tests not passing yet. I have to investigate the cause ...

@amartinhuertas
Copy link
Member Author

amartinhuertas commented Apr 12, 2021

Tests not passing yet. I have to investigate the cause ...

Ok. I figure out the cause. The tests which are failing are those corresponding to (oriented) simplices meshes.

The previous code was using the signed determinant when computing the face RT DoFs. See https://github.com/gridap/Gridap.jl/pull/577/files#diff-ba506509a89b9ff975762ef98f39e3dd9a7d0ef2b6fdb3b962576fd48abdd161L106

Thus, I guess that I have to also use cell_orientation for the evaluation of DoFs to fix this issue (this would match with what it seems they are doing in UFL/Dolphin). However, I still do not understand why. When I do hand-written calculations I do not see why I should change the sign of the DoF. (most probably I am missing something).

@fverdugo
Copy link
Member

Hi @amartinhuertas !

What is the difference between the cell_orientation array you mention and the result of get_cell_permutations(top::GridTopology,d::Integer) ? We use the output of the latter in several places to handle non-oriented meshes.

@fverdugo
Copy link
Member

fverdugo commented Apr 12, 2021

BTW, a better name for get_cell_permutations would be get_face_permutation to be consistent with the names of other routines.

@amartinhuertas
Copy link
Member Author

amartinhuertas commented Apr 12, 2021

What is the difference between the cell_orientation array you mention and the result of get_cell_permutations(top::GridTopology,d::Integer) ? We use the output of the latter in several places to handle non-oriented meshes.

I don't know. Can you elaborate a little bit more the output of get_cell_permutations(top::GridTopology,d::Integer) ?

I guess that you want to know if cell_orientation can be computed out of the output of get_cell_permutations, is that right?

cell_orientation is a cell-wise array, with as many entries as cells. As far as I understand, this array returns true for a given cell if the geometrical mapping from the reference cell to the physical cell involves a reflection transformation. For example, in 2D, the y and x axis of the coordinate system of the cell are flipped in physical space in case of a reflection transformation. In such a case, the determinant of the Jacobian is negative. However, in the case of manifolds, the pseudodeterminant does not carry a sign, thus we need another means to determine whether the coordinate system is flipped or not. This is the goal of the cell_orientation array. Does this make sense?

@fverdugo
Copy link
Member

fverdugo commented Apr 12, 2021

OK, at the end, it is just an array of bools, right? If so, I would name it cell_is_oriented (minor thing).

If these "cells" are "facets" of a volume mesh, I guess you can compute the array from the output of get_cell_permutataions as

topo = get_grid_topology(model)
D = num_cell_dims(model)
facet_is_oriented = lazy_map(i->i==1,get_cell_permutations(topo),D-1)

I.e., permutation index == 1 is for oriented facets.

Idem if the "cells" are the edges of a surface mesh.

@amartinhuertas
Copy link
Member Author

OK, at the end, it is just an array of bools, right? If so, I would name it cell_is_oriented (minor thing).

Yes. To be more precise, it is a cell-wise array of boolean Fields. I tried to have an array of bools instead, but not sure how I can achieve this given the structure of the code.

If these "cells" are "facets" of a volume mesh, I guess you can compute the array from the output of get_cell_permutataions as

These "cells" are the cells of a volume mesh or surface mesh embedded in a higher dim space, not facets of a volume mesh. In any case, I do not see why what you say is correct, (perhaps because I do not fully know the rationale underlying get_cell_permutations). My concern here is that you may have a value distinct from 1 for cell permutation without having a reflection transformation involved, right?

@fverdugo
Copy link
Member

Perhaps we can discuss this face-to-face!

  1. CellFE constructor is now invoked with cell_orientation.
  2. cell_orientation is now taken into account when
     computing the RT dofs.
@amartinhuertas
Copy link
Member Author

Perhaps we can discuss this face-to-face!

Sure! I guess that this would be very helpful. When are you available? Perhaps @santiagobadia should join as well ... he is aware of this development as well.

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
@amartinhuertas
Copy link
Member Author

The discussion above regarding cell_orientation + general manifolds no longer applies. I have undone all changes related to the cell_orientation array in this branch. Reason: see commit 276c410 description.

This PR depends on #579 and thus must be processed after it.

… github.com:gridap/Gridap.jl into misc_changes_rt_fes_on_general_manifolds

Conflicts:
	src/ReferenceFEs/RaviartThomasRefFEs.jl
	src/ReferenceFEs/ReferenceFEInterfaces.jl
@codecov-commenter
Copy link

codecov-commenter commented May 6, 2021

Codecov Report

Attention: Patch coverage is 92.30769% with 1 line in your changes missing coverage. Please review.

Project coverage is 87.66%. Comparing base (84fcf88) to head (81c62a4).
Report is 1825 commits behind head on master.

Files with missing lines Patch % Lines
src/ReferenceFEs/RaviartThomasRefFEs.jl 85.71% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #577      +/-   ##
==========================================
+ Coverage   87.64%   87.66%   +0.01%     
==========================================
  Files         132      132              
  Lines       13698    13701       +3     
==========================================
+ Hits        12006    12011       +5     
+ Misses       1692     1690       -2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.


🚨 Try these New Features:

@amartinhuertas
Copy link
Member Author

Hi @fverdugo!

After PR #577, now we can proceed with this PR (i.e., #579).

The changes are minor, essentially, replace det by meas, inv by pinvJt, and https://github.com/gridap/Gridap.jl/pull/577/files#diff-ea5ed427dd3dd1fad920c1f2014e8950d5866705bde29f7278c6d93390ae72c1R140.

These changes are required in order to be able to pass the test that I have introduced here:

https://github.com/gridap/Gridap.jl/pull/577/files#diff-37219e6468e24d7c4b29903834b894aa426a861aa49ef64122068cbd47ceb5daR69

@amartinhuertas
Copy link
Member Author

@santiagobadia ... @fverdugo is happy with the changes introduced in this PR. Are u ok as well? If yes, I think we can accept the PR.

Copy link
Member

@santiagobadia santiagobadia left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All good

@santiagobadia santiagobadia merged commit 8cb22e2 into master May 7, 2021
@fverdugo fverdugo deleted the misc_changes_rt_fes_on_general_manifolds branch May 7, 2021 10:24
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants