GridapDistributed.jl/src/DivAndCurlConformingFESpaces.jl at 2ac5643761478f1e225cb5f2a45f9bac2f3268ab · gridap/GridapDistributed.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
const RTorNedelec = Union{RaviartThomas,Nedelec}

function FESpaces.FESpace(model::DistributedDiscreteModel,
                          reffe::Tuple{<:RTorNedelec,Any,Any};
                          conformity=nothing,
                          split_own_and_ghost=false,
                          constraint=nothing,
                          kwargs...)

  cell_reffes = map(local_views(model)) do m
    basis,reffe_args,reffe_kwargs = reffe
    cell_reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
  end
  _common_fe_space_constructor(model,cell_reffes;conformity,split_own_and_ghost,kwargs...)
end

function FESpaces.FESpace(model::DistributedDiscreteModel,
                 reffe::GenericRefFE{<:RTorNedelec};
                 conformity=nothing,
                 split_own_and_ghost=false,
                 constraint=nothing,
                 kwargs...)
  cell_reffes = map(local_views(model)) do m
    Fill(reffe,num_cells(m))
  end
  _common_fe_space_constructor(model,cell_reffes;conformity,split_own_and_ghost,kwargs...)
end

function _setup_dmodel_and_dtrian(_trian::DistributedTriangulation)
  trian = add_ghost_cells(_trian)
  models = map(local_views(trian)) do t
    get_active_model(t)
  end
  GenericDistributedDiscreteModel(models, generate_cell_gids(trian)), trian
end

function FESpaces.FESpace(
  _trian::DistributedTriangulation,
  reffe::Tuple{<:RTorNedelec,Any,Any};
  conformity=nothing,
  split_own_and_ghost=false,
  constraint=nothing,kwargs...
)
  dmodel, dtrian = _setup_dmodel_and_dtrian(_trian)
  cell_reffes = map(local_views(dmodel)) do m
    basis,reffe_args,reffe_kwargs = reffe
    cell_reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
  end
  _common_fe_space_constructor(dmodel,cell_reffes,dtrian;conformity,split_own_and_ghost,kwargs...)
end

function FESpaces.FESpace(_trian::DistributedTriangulation,
                 reffe::GenericRefFE{<:RTorNedelec};
                 conformity=nothing,
                 split_own_and_ghost=false,
                 constraint=nothing,
                 kwargs...)
  dmodel, dtrian = _setup_dmodel_and_dtrian(_trian)
  cell_reffes = map(local_views(dmodel)) do m
     Fill(reffe,num_cells(m))
  end
  _common_fe_space_constructor(dmodel,cell_reffes,dtrian;conformity,split_own_and_ghost,kwargs...)
end

function _common_fe_space_constructor(model,cell_reffes,trian;conformity,split_own_and_ghost,kwargs...)
  sign_flips=_generate_sign_flips(model,cell_reffes)
  spaces = map(local_views(model),sign_flips,cell_reffes,local_views(trian)) do m,sign_flip,cell_reffe,trian
     conf = Conformity(testitem(cell_reffe),conformity)
     cell_fe = CellFE(m,cell_reffe,conf,sign_flip)
     FESpace(m, cell_fe; trian=trian, kwargs...)
  end
  gids = generate_gids(model,spaces)
  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
end

function _common_fe_space_constructor(model,cell_reffes;conformity,split_own_and_ghost,kwargs...)
  sign_flips=_generate_sign_flips(model,cell_reffes)
  spaces = map(local_views(model),sign_flips,cell_reffes) do m,sign_flip,cell_reffe
     conf = Conformity(testitem(cell_reffe),conformity)
     cell_fe = CellFE(m,cell_reffe,conf,sign_flip)
     FESpace(m, cell_fe; kwargs...)
  end
  gids = generate_gids(model,spaces)
  trian = DistributedTriangulation(map(get_triangulation,spaces),model)
  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
end

function _generate_sign_flips(model, cell_gids, cell_reffes)
sign_flips = map(local_views(model),partition(cell_gids),cell_reffes) do m, p, cell_reffe
    D = num_cell_dims(model)

    gtopo = get_grid_topology(m)

    # Extract composition among cells and facets
    cell_wise_facets_ids = get_faces(gtopo, D, D - 1)
    cache_cell_wise_facets_ids = array_cache(cell_wise_facets_ids)

    # Extract cells around facets
    cells_around_facets = get_faces(gtopo, D - 1, D)
    cache_cells_around_facets = array_cache(cells_around_facets)

    ncells = num_cells(m)
    ptrs = Vector{Int32}(undef,ncells+1)
    for cell in 1:ncells
      reffe=cell_reffe[cell]
      ptrs[cell+1] = num_dofs(reffe)
    end
    PArrays.length_to_ptrs!(ptrs)
    ndata = ptrs[end]-1
    data  = Vector{Bool}(undef,ndata)
    data .= false

    loc_to_glo=local_to_global(p)
    for cell in own_to_local(p)
      sign_flip = view(data,ptrs[cell]:ptrs[cell+1]-1)
      reffe = cell_reffe[cell]
      D = num_dims(reffe)
      face_own_dofs   = get_face_own_dofs(reffe)
      facet_lid       = get_offsets(get_polytope(reffe))[D] + 1
      cell_facets_ids = getindex!(cache_cell_wise_facets_ids,
                                  cell_wise_facets_ids,
                                  cell)
      for facet_gid in cell_facets_ids
          facet_cells_around = getindex!(cache_cells_around_facets,
                cells_around_facets,
                facet_gid)
          is_slave=false
          if length(facet_cells_around) > 1
            mx=maximum(loc_to_glo[facet_cells_around])
            is_slave = (loc_to_glo[cell] == mx)
          end
          if is_slave
              for dof in face_own_dofs[facet_lid]
                  sign_flip[dof] = true
              end
          end
          facet_lid = facet_lid + 1
      end
    end
    JaggedArray(data,ptrs)
  end
  cache = fetch_vector_ghost_values_cache(sign_flips,partition(cell_gids))
  fetch_vector_ghost_values!(sign_flips,cache) |> wait
  sign_flips
end


function _generate_sign_flips(model::DistributedDiscreteModel,
                             cell_reffes::AbstractVector{<:AbstractVector{<:GenericRefFE{<:RaviartThomas}}})
  cell_gids  = get_cell_gids(model)
  _generate_sign_flips(model,cell_gids,cell_reffes)
end

function _generate_sign_flips(model::DistributedDiscreteModel,
                              cell_reffes::AbstractVector{<:AbstractVector{<:GenericRefFE{<:Nedelec}}})
  Dc=num_cell_dims(model)
  is_simp = false
  map(cell_reffes) do cell_reffes
    cell_reffe = cell_reffes[1]
    polytope = get_polytope(cell_reffe)
    if (Gridap.Geometry.is_simplex(polytope))
        is_simp=true
    end
  end
  if (Dc==2 || is_simp)
    map(local_views(model), cell_reffes) do model, cell_reffes
        Gridap.FESpaces._no_sign_flip(model, cell_reffes)
    end
  else
    @assert Dc==3
    cell_gids  = get_cell_gids(model)
    _generate_sign_flips(model,cell_gids,cell_reffes)
  end
end