Higher order triangulation of sphere

[alecontri]

Hi,

I am trying to solve a Vector-Laplace equation on the sphere using Gridap. I wanted to test the influence of the geometry approximation, reason why I would like to use higher order geometries. I tried to look for possible ways to implement it but still haven't found a feasible way to do it. In particular:

• Generating a surface mesh with gmsh and trying to import the .msh file with the function "GmshDiscreteModel" from GridapGmsh throws the error ERROR: Only one element type per dimension allowed for the moment. Dimension 3 has 0 different element types. I attach a .txt file containing the text of the .geo file I am using sphere.txt
• Generating a volume mesh with gmsh, importing it using "GmshDiscreteModel" and then use only the boundary works, but it seems like it is not possible to import higher order .msh files and the error ERROR: For the moment only for first-order elements is thrown (Also, in this case how can I access normals?)
• Trying to use the ready made functions in "GridapGeosciences", precisely "https://github.com/gridapapps/GridapGeosciences.jl/blob/master/src/CubedSphereDiscreteModels.jl" and "https://github.com/gridapapps/GridapGeosciences.jl/blob/master/src/CubedSphereTriangulations.jl" seems to successfully complete the Triangulation, but when using "get_normal_vector" it throws ERROR: MethodError: get_normal_vector(::Gridap.Geometry.BodyFittedTriangulation{2, 3, Gridap.Geometry.UnstructuredDiscreteModel{2, 3, Float64, Gridap.Geometry.Oriented}, Gridap.Geometry.UnstructuredGrid{2, 3, Float64, Gridap.Geometry.Oriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}}) is ambiguous. I attach a minimal example using the above mentioned files: test.txt.

I was wondering if you know possible issues with the code or a simpler way to implement it.

I am using Julia v1.8.5, Gridap v0.17.17, GridapGmsh v0.6 and FillArrays v0.12.8.

Thanks a lot for the help

[amartinhuertas]

Hi @alecontri,

I would clearly go for the third option, i.e., the cubed sphere grid of the sphere (https://www.sciencedirect.com/science/article/abs/pii/S0021999196900479). This kind of mesh is pretty standard in numerical weather prediction cores, such as, e.g., LFric from the UkMet office. GridapGeosciences provides two possible geometrical mappings to represent the geometry (i.e., the mapping among 2D quadrilaters in reference space and quadrilaterals in 3D ambient space): (1) polynomial-based (i.e., FE-based) mapping of arbitrary order. (2) analytical.

I will take a look at the error you report with the get_normal_vector function. In any case:
(1) why do you need the normal? The following driver solves the Laplace-Beltrami surface PDE on the sphere (i.e. scalar-valued Laplacian), and I did not require the normal (https://github.com/gridapapps/GridapGeosciences.jl/blob/master/test/LaplaceBeltramiCubedSphereTests.jl).
(2) please note there are actually two possible normals: (a) the unit outward normal pointing in the radial direction. (b) the unit outward normal to each edge of the cube sphere cells, which belongs to the tangent space of the discrete manifold; this one is, e.g., the one used to evaluate the DoFs of an analytical function which you want to interpolate into the Raviart-Thomas vector-valued FE space.

[alecontri]

Hi @amartinhuertas,

thanks a lot for the prompt answer! I need the normal since I would like to solve the Vector Laplacian using surface finite elements as described for example in https://arxiv.org/abs/1610.06747 and I have to penalise the outward normal pointing component (in the radial direction) since tangentiality is not guarateed as in the Raviart-Thomas case.