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

Reporting a series of errors with void parts + tentative solutions #808

Closed
amartinhuertas opened this issue Jun 27, 2022 · 1 comment
Closed
Assignees

Comments

@amartinhuertas
Copy link
Member

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

@JordiManyer
Copy link
Member

I think this has been resolved.

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

No branches or pull requests

3 participants