gridap/Gridap.jl

Reporting a series of errors with void parts + tentative solutions

amartinhuertas opened this issue · 1 comments

Hi @fverdugo,

while exploring in a parallel context several scenarios in which the portion of a given processor is void, I have found 3 different errors, which I report below as MWEs along with possible solutions. (Note there is also an attached file needed to reproduce the first error).

Can you take a look at these, and let me know what do you think?

Once we agree on the approach to solution, I will send a PR with the corresponding fixes.

Thanks!
@amartinhuertas

using Gridap
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Io
using Gridap.FESpaces
using Gridap.CellData

u(x)=x[1]

# Fixes Err1 (see below)
function Gridap.Io.from_dict(::Type{UnstructuredGrid},dict::Dict{Symbol,Any})
  x = collect1d(dict[:node_coordinates])
  T = eltype(x)
  Dp = dict[:Dp]
  if (length(x)!=0)
    # Had to remove type annotation from this line (see corresponding line in Gridap)!!!
    # There seems to be a BUG in the compiler (Julia v1.7.2)
    node_coordinates = reinterpret(Point{Dp,T},x)
  else
    T1=Point{Dp,Float64}
    node_coordinates=Vector{T1}(undef,0)
    println(eltype(node_coordinates))
  end
  cell_node_ids = from_dict(Table{Int32,Vector{Int32},Vector{Int32}},dict[:cell_node_ids])
  reffes = [from_dict(LagrangianRefFE,reffe) for reffe in dict[:reffes]]
  cell_type::Vector{Int8} = dict[:cell_type]
  O::Bool = dict[:orientation]
  UnstructuredGrid(
    node_coordinates,
    cell_node_ids,
    reffes,
    cell_type,
    O ? Oriented() : NonOriented())
end

# Fixes Err2 (see below)
function Base.map(::typeof(testitem),
                  a::Tuple{<:AbstractVector{<:AbstractVector{<:VectorValue}},<:AbstractVector{<:Gridap.Fields.LinearCombinationFieldVector}})
  a2=testitem(a[2])
  a1=Vector{eltype(eltype(a[1]))}(undef,size(a2,1))
  a1.=zero(testitem(a1))
  (a1,a2)
end

# Fixes Err3 (see below)
function Gridap.Geometry.is_change_possible(
     strian::Gridap.Geometry.Triangulation,
     ttrian::Gridap.Geometry.Triangulation)
  if strian === ttrian || num_cells(strian)==num_cells(ttrian)==0
    return true
  end
  Gridap.Helpers.@check get_background_model(strian) === get_background_model(ttrian) "Triangulations do not point to the same background discrete model!"
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  is_change_possible(sglue,tglue) # Fails here
end

# Fixes Err3 (see below)
function Gridap.CellData.change_domain(a::CellField,
                                       ::ReferenceDomain,
                                       ttrian::Triangulation,
                                       ::ReferenceDomain)
  msg = """\n
  We cannot move the given CellField to the reference domain of the requested triangulation.
  Make sure that the given triangulation is either the same as the triangulation on which the
  CellField is defined, or that the latter triangulation is the background of the former.
  """
  strian = get_triangulation(a)
  if strian === ttrian || num_cells(strian)==num_cells(ttrian)==0
    return a
  end
  @assert Gridap.Geometry.is_change_possible(strian,ttrian) msg
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  change_domain_ref_ref(a,ttrian,sglue,tglue)
end

model=Gridap.Io.from_json_file(DiscreteModel,"model_part_2.json") # Fails here (Err1)

# FE Spaces
order=1
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe,dirichlet_tags="boundary")
U = TrialFESpace(u,V) # Fails here (Err2)

s     = get_fe_dof_basis(V)
trian = get_triangulation(s)
f     = CellField(u,trian,DomainStyle(s))
b = change_domain(f,s.domain_style)
s(f) # Fails here (Err2 at a lower level)

sd=get_data(s)
bd=get_data(b)
sdi=testitem(sd)
sbd=testitem(bd)
evaluate(sdi,sbd) # Fails here (Err2 at a lower level)

uh = interpolate(u,U)
Ω  = view(get_triangulation(U),Int[])
dΩ = Measure(Ω,2)
dc=(uh)dΩ # Fails here (Err3)

trian_f=get_triangulation(uh)
trian_x=.quad.trian
Gridap.CellData.is_change_possible(trian_f,trian_x) # Fails here (Err3)

model_part_2.zip

I think this has been resolved.