Future developments about gridap.jl

[santiagobadia]

@fverdugo

As suggested, I have added two new methods to Polytope

One is function generate_admissible_permutations(p::Polytope), which provides the list of possible enumerations of a polytope. An admissible enumeration is uniquely determined by the first D+1 nodes.

I have also included a new method

function _dimfrom_fs_dimto_fs(p::Polytope, dim_from::Int, dim_to::Int)

which provides the n-faces of dim_to of all n-faces of dim_from in a given polytope. Example of usage:

p = Polytope(1,1,1)
nf_vs = Gridap.Polytopes._dimfrom_fs_dimto_fs(p,2,0)
@test length(nf_vs) == 6

[fverdugo]

Other things also worth considering:

• Automatic code style revision (i.e. look for trailing white space, tabs, max line width violations, etc. ). This can be integrated in the test suite.
• Meaningful error messages

[santiagobadia]

Some missing things:

• Use abstract interface getters in FESpaces module instead of access to concrete data
• Create test FE space from FE Space (w/ homogeneous Dirichlet data)
• Create trial FE space from FE Space (w/ non-homogeneous Dirichlet data)
• Interpolators for FE Space without Dirichlet data and the FE Spaces above with Dirichlet data, where the boundary values are respected in the first case but not in the other ones
• Create multi-field Assembler with an array of FESpaces
• Integrate a multi-field bilinear form with arrays of FEFunctions, etc
• Create const for FEFunction in CellMaps, and fixed/free DOF values getters
• Define FEOperator, which provides the residual for an array of FEFunctions (Vararg type)
• Implement an evaluate and solve method for FEOperator
• Consider problems with physical values that are also FEFunctions
• Use automatic differentiation to obtain the system matrix from the residual value (cell-wise or global?)
• Put in FESpace an array of PhysicalTags that define the Dirichlet boundary, waiting for @fverdugo work
• Determine a Dirichlet data interpolator (for a FEFunction) in which for every PhysicalTag that is in the FESpace Dirichlet boundary array, we provide a vector of analytical fields
• Put an assert when a function to be interpolated provides an Int
• Extract the fixed dofs info from PhysicalTag in the FESpace constructor
• Create strucs NonConformingFESpace with ConstrainingDOFs and ConstrainingValues arrays
• Create Constrained FESpace that composes an abstract FESpace plus ConstrainingDOFs and ConstrainingValues arrays

After @fverdugo re-factorization of CellMaps

• Fix the problems with evaluate of cell vectors/matrices
• Create a MultiValuedCellValue that generalizes ConstantCellValue, and use it in FESpace, etc

Other misc (longer-term) issues:

• Make all tet grids oriented at the very beginning
• Implement change of basis for the DOF equivalence class for non-oriented meshes (e.g., unstructured hex meshes)

[fverdugo]

• Put in FESpace an array of PhysicalTags that define the Dirichlet boundary, waiting for @fverdugo work

Done!

[fverdugo]

After @fverdugo re-factorization of CellMaps

Refactoring of CellValues and CellMaps will take place in branch cell_values_and_cell_maps_refactoring

[santiagobadia]

@fverdugo : I would like to create a wrapper of FESpace that includes the Dirichlet Data. However, I don't want to replicate all the methods already defined for a concrete FESpace in the wrapper. I would like to have something like this:

struct FESpaceWithDirichletData{D,Z,T,E,V<:FESpace{D,Z,T,E}} <: V
  fesp::V
  dir_data::Vector{Float64}
end

This way FESpace + Dirichlet data would always sub-type the original FESpace concretization, and it would not require any method implementation for FESpaceWithDirichletData but dir_data related things.

Unsurprisingly, it does not work, wrong subtyping. I guess it cannot check what V is and whether it allows sub-typing, etc (?)

Any other solution?

[santiagobadia]

Probably, I can just do:

for op in (all FESpace interface methods)
  @eval begin
      $op(this::FESpaceWithDirichletData) = $op(this.fesp)
    end
end

since all current methods only have one FESpace as input arguments...

[santiagobadia]

It works!

We could probably use it to automatically generate getters for all structs by iterating over their attributes

It is not going to work for more complicated interfaces

[fverdugo]

In julia, it is only allowed to subtype abstract types, and surely there is a good reason.
In consequence, I would not try to emulate extension of concrete types...

Any way, even if it would be allowed to extend concrete types, I would not do it here since we want to extend with Dirichlet data ANY FE space (and not only a concrete implementation)...

I would do a concrete struct like this

struct FESpaceWithDirichletData{F<:AbstractFESpace} <: AbstractFESpace
  fespace::F
  diridata::TypeThatHoldsDiriData
end

Then FESpaceWithDirichletData will delegate all queries that are required to the internal fespace

[santiagobadia]

I finally did

struct FESpaceWithDirichletData{D,Z,T,E,V<:FESpace{D,Z,T,E}} <: FESpace{D,Z,T,E}
	fesp::V
	dir_data::Vector{Float64}
end

for op in (:reffes, :triangulation, :gridgraph, :nf_eqclass, :cell_eqclass,
	:num_free_dofs, :num_fixed_dofs)
	@eval begin
		$op(this::FESpaceWithDirichletData) = $op(this.fesp)
	end
end

My only concern was re-implementing all the getters, but it is not so bad

[fverdugo]

Yes! this is what I meant.

The way you avoid code repetition with the @eval is totally fine.

[santiagobadia]

@fverdugo

I want to create a wrapper for FESpaces that adds dirichlet data to it, as we have commented many times. However, if my FESpace is of type X, I would like that the FESpaceWithDirichletData could be re-used in the functions defined for X type, i.e., the interpolation function for a ConformingFESpace should be also used for the FESpaceWithDirichletData that is generated with a ConformingFESpace.

I was thinking on something like this

abstract type str{T} end

struct str1{T} <: str{T}
	a::T
end

struct str2{T,S<:str{T}}
	a::T
	s::S
end

strX{T} = Union{str1{T},str2{T,str1{T}}}

What do you think? str should be replace by FESpace, etc...

[fverdugo]

This is the typical situation where "traits" are useful.

For instance:

# we define here the trait
abstract type ConformityStyle end
struct ConformingStyle <: ConformityStyle end
struct NonConformingStyle <: ConformityStyle end

abstract type AbstractFESpace end

# Each FE space has to define which is its conformity style
ConformityStyle(::AbstractFESpace)::ConformityStyle = @abstractmethod 

# insted of an abstract method one can also use non-conforming by default:
ConformityStyle(::AbstractFESpace) = NonConformingStyle()

then imagine that there is a query (e.g., enum_dofs) that can be optimized depending on the conformity, one can implement two variants:

function enum_dofs(f::AbstractFESpace,::NonConformingStyle)
# Implementation for **ANY** FEspace that is NonConforming
...
end

function enum_dofs(f::AbstractFESpace,::ConformingStyle)
# Implementation for **ANY** FEspace that is Conforming
...
end

# These functions are used as follows

fespace = MyFESpaceConcreteType()

enum_dofs( fespace, ConformityStyle(fespace) )

This is very usefull also if we want to dispatch on multiple traits, e.g., one could also have other traints like OrderStyle, (FirstOrderStyle, HighOrderStyle), HomogeneityStyle (SingleRefFEStyle, MultipleRefFEStyle), which avoids ending with types like ConfonformgAndFirstOrderFESpace, ConformingAndHighOrderFESpace ...

Julia Base is full of this pattern. e.g. IndexStyle for arrays... If you think it is useful in your situation, use use traits.

[fverdugo]

I think that it makes sense to have this pattern in FESpace, since in the future we want to optimize for different situations. The other good thing is that you do not have to anticipate for which traits you want to optimize (adding a new trait in the future is straight forward)

[santiagobadia]

At the end, it is just to add a new element in the type hierarchy.
I can do it this way, it is simple.

[fverdugo]

Actually, traits are also useful since you cannot extend from concrete types. They allow to "inherit" the optimizations when you wrap a type.

[fverdugo]

Also useful since there is no multiple inheritance

[santiagobadia]

I will look at this now. I have read some discussions on it.
At the end of the day, it defines a subset of types independent of the type hierarchy.
It is pretty cool.

[santiagobadia]

I am not sure traits work here, since it seems that one type must have one trait, and in my case, I want it to be different, depending on the values of the type. In your previous example, I think that it should be

ConformityStyle(::Type{AbstractFESpace})::ConformityStyle = @abstractmethod 

instead of

ConformityStyle(::AbstractFESpace)::ConformityStyle = @abstractmethod 

So, it seems it does not solve my problem... I will do more research

[santiagobadia]

Ok, using templatization I can fix it

[santiagobadia]

abstract type mystyle end
struct style1 <: mystyle end
struct style2 <: mystyle end

abstract type str{T} end
mystyle(::str)::mystyle = error("Not defined")

struct str1{T} <: str{T}
	a::T
end
mystyle(::Type{str1{T}}) where {T}= style1

struct str2{T} <: str{T}
	a::T
end
mystyle(::Type{str2{T}}) where {T} = style2

struct strc{T,S<:str{T}}
	a::T
	s::S
end
mystyle(::Type{strc{T,S}}) where {T,S<:str1{T}} = style1
mystyle(::Type{strc{T,S}}) where {T,S<:str2{T}} = style2

s1 = str1(10)
s2 = str2(10)
s1c = strc(20,s1)
s2c = strc(20,s2)



mystyle(typeof(s1)) # style1
mystyle(typeof(s2)) # style2
mystyle(typeof(s1c)) #style1
mystyle(typeof(s2c)) #style2

[fverdugo]

Yes, the style getter applyes in the type instead of the value.

If you want to get the style of the wrapped object type, it can also be done. Actually we are already doing it in some places. Eg:
https://github.com/lssc-team/Numa.jl/blob/83bfd3f532430437548494670763f3a2e8950c8f/src/CellValues/Wrappers.jl#L37

[santiagobadia]

Yes, it is an obvious wrapper for the type-based function.
What I meant before was that I do not think you want to implement something like

mystyle(this::strc) = mystyle(strc.s)

because we would not be able to implement the type-based function (it is not well defined).

[santiagobadia]

OK, I understand now what you ment. It gives the same result, but it is better.

[santiagobadia]

Implemented... I am using a wrapper function to not to show it to users.

Any syntatic sugar in Julia? I will take a look...

function interpolate(fun::Function, fesp::ConformingFESpaces{D}) where {D}
	@assert MeshConformity(fesp) == ConformingMesh() # Not implemented
	_interpolate(fun,fesp, MeshConformity(fesp))
end

function _interpolate(fun::Function, fesp::ConformingFESpaces{D}, ::ConformingMesh) where {D}
...

[fverdugo]

Target driver single field

model = GMeshModel("model.msh")
g1(x) = 2*x[1]
g2(x) = 3*x[2]
f1(x) = 5*x[3]
h1(x,u) = 3*x[3]+0.5*u
ufun(x) = 3*x[2]
gradufun(x) = VectorValue(0.0,2.0,0.0)
gradient(::typeof(ufun)) = gradufun

trian = Triangulation(model)
quad = CellQuadrature(trian,order=2)

fespace = FESpace(model,order=1,diritags=[3,34])
V = TestSpace(fespace)
U = TrialSpace(fespace,[g1,g2])

f = CellField(f1,trian)
nu(u) = CellField(h1,trian,u)

dr(u,v,du) = inner(gradient(v), nu(u)*gradient(du) )

r(u,v) = inner(v,f) - dr(u,v,u)

A = FEOperator(r,dr,V,U,quad)

uh = solve(A,tol=10e-3)

uexa = CellField(ufun,trian)

e = uexa - uh

a(u,u) = dr(u,u,u)

err = integrate(a(e,e),trian,quad)
err = sqrt( sum(err) )
@test err < 10e-10

writevtk(trian,"fe_solution",cellfields=["uh"=>uh])