Develop a new NodesArray for general polytopes

[fverdugo]

@santiagobadia

It seems that there is a problem in NodesArray for simplices polytopes:

using Gridap
const t = TET_AXIS
polytope = Polytope((t,t)) # This is a triangle
orders=[1,1]
na = NodesArray(polytope,orders)
length(na.coordinates) == 4 # Wrong? Should be 3, right?

The length of na.coordinates seems to be wrong. Could you confirm that this is a bug?

[santiagobadia]

More than a typo, it is a missing development. The nodearray is for n-cubes only.
It is an essential development to be completed asap.
I have been trying to keep the polytope in good shape but I have not touched the NodesArray for ages.

[fverdugo]

Ok! Then we need a not implemented error.
I have added it in line

Gridap.jl/src/RefFEs/Polytopes.jl

Line 105 in 785eb61

@notimplementedif any([ t != HEX_AXIS for t in polytope.extrusion])

[fverdugo]

I remove also the bug flag since now it is only a functionality to be added

[fverdugo]

I'll push into master some changes that allow one to construct a grid from a polytope. From the grid one can build a triangulation as always, and from this a CellGeomap. With this CellGeomap one can map nodes from a low dim polytope to a higher dim one:

t = HEX_AXIS
polytope = Polytope((t,t,t))
grid = Grid(polytope,2) # A grid formed by the faces of the hexahedron (requires a NodesArray of order 1, which should be easy to construct in general)
trian = Triangulation(grid)
quad = CellQuadrature(trian,order=2)
q = coordinates(quad) # if you replace this with the nodes on the reference face, you will be almost done for the NodesArray
phi = CellGeomap(trian)
x = evaluate(phi,q)
writevtk(grid,"grid")
writevtk(x,"x")

With this, implementing a high order NodesArray for a general polytope should be easy.

[santiagobadia]

@fverdugo I have created a new method generate_interior_nodes in Polytope that given an order returns an array of interior nodes for an n-face. I have stored the nodes as Point{Int} with coordinates that go from 1 to order-1

I think that with this method we can close the issue, since I could readily generate all the nodes easily.

The boundary nodes will be computed with lower dimension n-faces using the previous method plus geomap. The geomap part is missing.

[fverdugo]

Great!

Now, we need another method generate_vertex_nodes (or whatever) that generates the coordinates of the vertices of the poltytope. Perhaps, this is already implemented with a different API.

Then, we can combine the generate_interior_nodes with generate_vertex_nodes using the CellGeomap as commented above to create a high order NodesArray. Here, one has to be careful since in order to construct the CellGeomap one needs to create a linear lagrangian ref fe. This ref fe cannot use the machinery of high order NodesArray for the degenerated case order = 1 since we would end up in a dependency loop. Perhaps, linear Lagrangian reffes and high order ones will need to be implemented separately for this reason. The linear ones using the generate_vertex_nodes and the high order ones using the high order NodesArray. I need to think more about this...

I would close the issue when everything is implemented. Makes sense, right?

[fverdugo]

I also forgot to mention that we also need to rescale the coordinates given by generate_interior_nodes to the interval [-1,1] or [0,1]. We take [0,1] for all polytopes by convenion

Just a question @santiagobadia , the argument p of generate_interior_nodes has to be a NFace ? I think that it would be easier for the user code if it is a Polytope instead of a NFace. Change NFace -> Polytope

[fverdugo]

• implement generate_interior_nodes returning Point{D,Float64} in [0,1]^D
• implement generate_vertex_nodes returning Point{D,Float64} in [0,1]^D
• change NFace -> Polytope in generate_interior_nodes
• Nodes array for the closed polytope
• Implement method to determine whether inwards or outwards normal for each facet

[santiagobadia]

I have added new methods in polytope, equidistant_interior_nodes_coordinates to get the coordinates of the nodes in the interior of a polytope, and vertices_coordinates to get the array of vertices.

Next, in the RefFE, for linear elements I have created a method that provides the old NodesArray info in _linear_lagrangian_nodes_polytope. For high order FEs, I have created _high_order_lagrangian_nodes_polytope, which goes through all n-faces, generate the interior nodes of the reference polytope of the n-face, and (using a linear RefFE for this polytope), maps it into the n-face using _map_high_order_lagrangian_nodes.

All these methods are being used in the RefFE constructor.

Now, we can simply write

p = Polytope(HEX_AXIS,TET_AXIS,TET_AXIS)
LagrangianRefFE{3,Float64}(p,3)
grid = Grid(p,2)
writevtk(grid,"grid")

to get a P3 FE on a tet.

I don't know why, but Grid print is not working. @fverdugo can you take a look. I have replaced the old NodesArray code with the new methods in Grid and seems to work. But I get an error when printing.

For the moment I do not have a NodesArray. It is quite simple to create one, but I am not sure it is needed. In fact, it was not in RefFE. Adding it should be straighforward. For the moment, I keep the old one because for anisotropic order the new version is not working yet, we consider same order in all dimensions. On the other hand, the generation of monomials is only working for n-cubes and n-tets. We should probably think about changes in the polynomial machinery to include other cases.

[santiagobadia]

@fverdugo when you take a look at the printing issue, we can probably close this issue. The remaining points are in #49

[santiagobadia]

In the way to eliminate NodesArrays I did some simple changes in GeometryWrappers.jl, but after them, some integration on faces tests do not work. Thus, the new version is in GeometryWrappersNew.jl and I have kept the old one with NodesArray. Is it probably related to the [-1,1] to [0,1] change?

[fverdugo]

@fverdugo when you take a look at the printing issue, we can probably close this issue. The remaining points are in #49

The implementation of the mask providing information about to where the normal vector points is also to be done, right? or it is already done?

[fverdugo]

In the way to eliminate NodesArrays I did some simple changes in GeometryWrappers.jl, but after them, some integration on faces tests do not work. Thus, the new version is in GeometryWrappersNew.jl and I have kept the old one with NodesArray. Is it probably related to the [-1,1] to [0,1] change?

Yes, we have to rescale the 1d quadrature if we have not done yet so.

[santiagobadia]

Done, also in old NodesArray

[santiagobadia]

With regard to the normals, I have created a method that provides the outward normals (checked for 2 to 4D and hex and tets).
The function is facet_normals in Polytopes.jl. Is that OK? The +1/-1 could also be computed, it is in fact part of the computation.

[fverdugo]

In fact, I only need (at least for the moment) the one that provides +1/-1. So, it would be nice to have this in a separate function. Is it possible?

[santiagobadia]

Ok, added this info, now facet_normals provide the normals and facet orientations (+1,-1)

D = 3
p = Polytope(1,1,1)
ns, f_os = Gridap.Polytopes.facet_normals(p)

[santiagobadia]

@fverdugo, Nice picture ;-)

I am not sure about the orientation I provide.

In a first approach I used the cross product for 3D, which naturally provides the orientation I think you need. However, the cross product only works in 3D and 7D. I think we needed something more general. Instead, I decided to compute the normal vector as the nullspace of a matrix with rows the vectors from all vertices but the first one to the first vertex of a facet. Next, I take the vector from a vertex of the polytope that does not belong to the facet to a vertex that belongs to it (the first one) in order for the normal vector to point outwards. This way, we have a clean and general way to determine the normal for any dimension and any polytope.

However, now the concept of orientation I provide does not have much sense (I think). In DG methods you can compute everything by transforming your outward normal in the reference space to the physical one. In many situations the normal is needed (think about DG methods for Darcy or Maxwell, etc). And it can also be used for any DG methods. No additional information is needed.

We can probably mathematically define the orientation you need for a general polytope and dimension and I can try to implement it. In the meantime, I think I will eliminate the orientation.

[fverdugo]

Hi @santiagobadia, I am sitting in the train returning from valencia.

I think, what we are looking for is the definition of outward normal in the sense of the divergence therorem.

Do you know any statement of this theorem in arbitrary dimensions?

This paper is for the stokes theorem in higher dims:

http://math.uchicago.edu/~may/REU2015/REUPapers/Kurila.pdf

[santiagobadia]

@fverdugo there is only one outwards normal, and it is the one I provide.

I think what you want is to stick a normal on a reference polytope for a face, e.g., the [1,0]x{0} segment for the faces of a square and n=[0,1]. Next map the reference face to any face and put 1 if the resulting normal points outwards and -1 if not. Is this what you need?

I think we can talk next week because with the orientation array (for gluing DOFs, GPs) and the outwards normal, I believe we should be able to implement any DG method.

We can probably talk to clarify it.

[fverdugo]

@santiagobadia I have fixed the conversion of Polytope to Grid.

Now the following code should work

p = Polytope(TET_AXIS,TET_AXIS,TET_AXIS) # Note that the first axis is TET_AXIS aswell, not HEX_AXIS as in your example 
grid = Grid(p,2)
writevtk(grid,"grid")

It is important to note that we need that all axis are TET_AXIS. This is because the conversion currently assumes that the faces of the resulting grid have the same extrusion code (see issue #50 ). For tets and hexs this is enough for the moment.

On the other hand, I have seen that you have rescaled all domains from [-1,1] to [0,1] as we wanted. This is an important change that needs to be reflected in NEWS.md (I have already included it).

I close the issue. For the discussion on the normals, I have created a new issue (see issue #51), where we can continue the discussion.