MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ptransfermap.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "psubmesh.hpp"
17 #include "ptransfermap.hpp"
18 #include "submesh_utils.hpp"
19 
20 using namespace mfem;
21 
23  const ParGridFunction &dst)
24 {
25  const ParFiniteElementSpace *parentfes = nullptr, *subfes1 = nullptr,
26  *subfes2 = nullptr;
27 
30  {
31  ParSubMesh* src_sm = static_cast<ParSubMesh*>(src.ParFESpace()->GetParMesh());
32  ParSubMesh* dst_sm = static_cast<ParSubMesh*>(dst.ParFESpace()->GetParMesh());
33 
34  // There is no immediate relation and both src and dst come from a
35  // SubMesh, check if they have an equivalent root parent.
36  if (SubMeshUtils::GetRootParent<ParSubMesh, ParMesh>(*src_sm) !=
37  SubMeshUtils::GetRootParent<ParSubMesh, ParMesh>(*dst_sm))
38  {
39  MFEM_ABORT("Can't find a relation between the two GridFunctions");
40  }
41 
43 
44  root_fes_ = new ParFiniteElementSpace(*src.ParFESpace(),
45  *const_cast<ParMesh *>(SubMeshUtils::GetRootParent<ParSubMesh, ParMesh>
46  (*src_sm)));
47  subfes1 = src.ParFESpace();
48  subfes2 = dst.ParFESpace();
49 
51  *root_fes_,
52  src_sm->GetFrom(),
53  src_sm->GetParentElementIDMap(),
54  sub1_to_parent_map_);
55 
57  *root_fes_,
58  dst_sm->GetFrom(),
59  dst_sm->GetParentElementIDMap(),
60  sub2_to_parent_map_);
61 
62  root_gc_ = &root_fes_->GroupComm();
63  CommunicateIndicesSet(sub1_to_parent_map_, root_fes_->GetVSize());
64 
65  z_.SetSize(root_fes_->GetVSize());
66  }
68  {
70  ParSubMesh* src_sm = static_cast<ParSubMesh*>(src.ParFESpace()->GetParMesh());
71  subfes1 = src.ParFESpace();
72  parentfes = dst.ParFESpace();
74  *parentfes,
75  src_sm->GetFrom(),
76  src_sm->GetParentElementIDMap(),
77  sub1_to_parent_map_);
78 
79  root_gc_ = &parentfes->GroupComm();
80  CommunicateIndicesSet(sub1_to_parent_map_, dst.Size());
81  }
83  {
85  ParSubMesh* dst_sm = static_cast<ParSubMesh*>(dst.ParFESpace()->GetParMesh());
86  subfes1 = dst.ParFESpace();
87  parentfes = src.ParFESpace();
89  *parentfes,
90  dst_sm->GetFrom(),
91  dst_sm->GetParentElementIDMap(),
92  sub1_to_parent_map_);
93  }
94  else
95  {
96  MFEM_ABORT("Trying to do a transfer between GridFunctions but none of them is defined on a SubMesh");
97  }
98 }
99 
101  ParGridFunction &dst) const
102 {
103  if (category_ == TransferCategory::ParentToSubMesh)
104  {
105  // dst = S1^T src
106  for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
107  {
108  dst(i) = src(sub1_to_parent_map_[i]);
109  }
110  }
111  else if (category_ == TransferCategory::SubMeshToParent)
112  {
113  // dst = G S1 src
114  // = G z
115  //
116  // G is identity if the partitioning matches
117 
118  for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
119  {
120  dst(sub1_to_parent_map_[i]) = src(i);
121  }
122 
123  CommunicateSharedVdofs(dst);
124  }
125  else if (category_ == TransferCategory::SubMeshToSubMesh)
126  {
127  // dst = S2^T G (S1 src (*) S2 dst)
128  //
129  // G is identity if the partitioning matches
130 
131  z_ = 0.0;
132 
133  for (int i = 0; i < sub2_to_parent_map_.Size(); i++)
134  {
135  z_(sub2_to_parent_map_[i]) = dst(i);
136  }
137 
138  for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
139  {
140  z_(sub1_to_parent_map_[i]) = src(i);
141  }
142 
143  CommunicateSharedVdofs(z_);
144 
145  for (int i = 0; i < sub2_to_parent_map_.Size(); i++)
146  {
147  dst(i) = z_(sub2_to_parent_map_[i]);
148  }
149  }
150  else
151  {
152  MFEM_ABORT("unknown TransferCategory: " << category_);
153  }
154 }
155 
156 void ParTransferMap::CommunicateIndicesSet(Array<int> &map, int dst_sz)
157 {
158  indices_set_local_.SetSize(dst_sz);
159  indices_set_local_ = 0;
160  for (int i = 0; i < map.Size(); i++)
161  {
162  indices_set_local_[map[i]] = 1;
163  }
164  indices_set_global_ = indices_set_local_;
165  root_gc_->Reduce(indices_set_global_, GroupCommunicator::Sum);
166  root_gc_->Bcast(indices_set_global_);
167 }
168 
169 void ParTransferMap::CommunicateSharedVdofs(Vector &f) const
170 {
171  // f is usually defined on the root vdofs
172 
173  const Table &group_ldof = root_gc_->GroupLDofTable();
174 
175  // Identify indices that were only set by other ranks and clear the dof.
176  for (int i = 0; i < group_ldof.Size_of_connections(); i++)
177  {
178  const int j = group_ldof.GetJ()[i];
179  if (indices_set_global_[j] != 0 && indices_set_local_[j] == 0)
180  {
181  f(j) = 0.0;
182  }
183  }
184 
185  // TODO: do the reduce only on dofs of interest
186  root_gc_->Reduce<double>(f, GroupCommunicator::Sum);
187 
188  // Indices that were set from this rank or other ranks have been summed up
189  // and therefore need to be "averaged". Note that this results in the exact
190  // value that is desired.
191  for (int i = 0; i < group_ldof.Size_of_connections(); i++)
192  {
193  const int j = group_ldof.GetJ()[i];
194  if (indices_set_global_[j] != 0)
195  {
196  f(j) /= indices_set_global_[j];
197  }
198  }
199 
200  // Indices for dofs that are shared between processors need to be divided by
201  // the whole group size that share this dof.
202  for (int gr = 1; gr < group_ldof.Size(); gr++)
203  {
204  for (int i = 0; i < group_ldof.RowSize(gr); i++)
205  {
206  const int j = group_ldof.GetRow(gr)[i];
207  if (indices_set_global_[j] == 0)
208  {
209  f(j) /= root_gc_->GetGroupTopology().GetGroupSize(gr);
210  }
211  }
212  }
213 
214  root_gc_->Bcast<double>(f);
215 }
216 
218 {
219  delete root_fes_;
220 }
221 
222 #endif // MFEM_USE_MPI
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
void BuildVdofToVdofMap(const FiniteElementSpace &subfes, const FiniteElementSpace &parentfes, const SubMesh::From &from, const Array< int > &parent_element_ids, Array< int > &vdof_to_vdof_map)
Build the vdof to vdof mapping between two FiniteElementSpace objects.
int * GetJ()
Definition: table.hpp:114
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
int GetGroupSize(int g) const
Get the number of processors in a group.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
Abstract parallel finite element space.
Definition: pfespace.hpp:28
SubMesh::From GetFrom()
Get the From indicator.
Definition: psubmesh.hpp:96
int Size_of_connections() const
Definition: table.hpp:98
Subdomain representation of a topological parent in another ParMesh.
Definition: psubmesh.hpp:51
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
static bool IsParSubMesh(const ParMesh *m)
Check if ParMesh m is a ParSubMesh.
Definition: psubmesh.hpp:170
void Transfer(const ParGridFunction &src, ParGridFunction &dst) const
Transfer the source ParGridFunction to the destination ParGridFunction.
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
const GroupTopology & GetGroupTopology()
Get a reference to the associated GroupTopology object.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
GroupCommunicator & GroupComm()
Return a reference to the internal GroupCommunicator (on VDofs)
Definition: pfespace.hpp:338
Table & GroupLDofTable()
Fill-in the returned Table reference to initialize the GroupCommunicator then call Finalize()...
ParTransferMap(const ParGridFunction &src, const ParGridFunction &dst)
Construct a new ParTransferMap object which transfers degrees of freedom from the source ParGridFunct...
Vector data type.
Definition: vector.hpp:60
int RowSize(int i) const
Definition: table.hpp:108
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
const Array< int > & GetParentElementIDMap() const
Get the parent element id map.
Definition: psubmesh.hpp:106
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113