GridapDistributed.jl/src/DivConformingFESpaces.jl at e02ed5179389668c67217f89fb3188c505869f0b · 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
function FESpaces.FESpace(model::DistributedDiscreteModel,
                          reffe::Tuple{RaviartThomas,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 FESpace(model::DistributedDiscreteModel,
                 reffe::GenericRefFE{RaviartThomas};
                 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 FESpace(
  _trian::DistributedTriangulation,
  reffe::Tuple{RaviartThomas,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 FESpace(_trian::DistributedTriangulation,
                 reffe::GenericRefFE{RaviartThomas};
                 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_reffes)
  cell_gids  = get_cell_gids(model)
  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)
            is_slave == false
          else
            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