MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ptransfermap.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
13
14#ifdef MFEM_USE_MPI
15
16#include "psubmesh.hpp"
17#include "ptransfermap.hpp"
18#include "submesh_utils.hpp"
19
20using 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(*src_sm) !=
38 {
39 MFEM_ABORT("Can't find a relation between the two GridFunctions");
40 }
41
43
44 {
45 ParMesh * parent_mesh =
46 const_cast<ParMesh *>(SubMeshUtils::GetRootParent(*src_sm));
47
48 int parent_dim = parent_mesh->Dimension();
49 int src_sm_dim = src_sm->Dimension();
50 int dst_sm_dim = dst_sm->Dimension();
51
52 bool root_fes_reset = false;
53 if (src_sm_dim == parent_dim - 1 && dst_sm_dim == parent_dim - 1)
54 {
55 const ParFiniteElementSpace *src_fes = src.ParFESpace();
56 const ParFiniteElementSpace *dst_fes = dst.ParFESpace();
57
58 const FiniteElementCollection *src_fec = src_fes->FEColl();
59 const FiniteElementCollection *dst_fec = dst_fes->FEColl();
60
61 const L2_FECollection *src_l2_fec =
62 dynamic_cast<const L2_FECollection*>(src_fec);
63 const L2_FECollection *dst_l2_fec =
64 dynamic_cast<const L2_FECollection*>(dst_fec);
65
66 if (src_l2_fec != NULL && dst_l2_fec != NULL)
67 {
68 // Source and destination are both lower dimension L2 spaces.
69 // Transfer them as the trace of an RT space if possible.
70
71 int src_mt = src_fec->GetMapType(src_sm_dim);
72 int dst_mt = dst_fec->GetMapType(dst_sm_dim);
73
74 int src_bt = src_l2_fec->GetBasisType();
75 int dst_bt = dst_l2_fec->GetBasisType();
76
77 int src_p = src_fec->GetOrder();
78 int dst_p = dst_fec->GetOrder();
79
80 if (src_mt == FiniteElement::INTEGRAL &&
81 dst_mt == FiniteElement::INTEGRAL &&
82 src_bt == BasisType::GaussLegendre &&
83 dst_bt == BasisType::GaussLegendre &&
84 src_p == dst_p)
85 {
86 // The subspaces are consistent with the trace of an RT space
87 root_fec_.reset(new RT_FECollection(src_p, parent_dim));
88 root_fes_.reset(new ParFiniteElementSpace(
89 const_cast<ParMesh *>(
91 root_fec_.get()));
92 root_fes_reset = true;
93 }
94 }
95 }
96
97 if (!root_fes_reset)
98 {
99 root_fes_.reset(new ParFiniteElementSpace(
100 *src.ParFESpace(),
101 const_cast<ParMesh *>(
102 SubMeshUtils::GetRootParent(*src_sm))));
103 }
104 }
105
106 subfes1 = src.ParFESpace();
107 subfes2 = dst.ParFESpace();
108
110 *root_fes_,
111 src_sm->GetFrom(),
112 src_sm->GetParentElementIDMap(),
113 sub1_to_parent_map_);
114
116 *root_fes_,
117 dst_sm->GetFrom(),
118 dst_sm->GetParentElementIDMap(),
119 sub2_to_parent_map_);
120
121 root_gc_ = &root_fes_->GroupComm();
122 CommunicateIndicesSet(sub1_to_parent_map_, root_fes_->GetVSize());
123
124 z_.SetSize(root_fes_->GetVSize());
125 }
127 {
129 ParSubMesh* src_sm = static_cast<ParSubMesh*>(src.ParFESpace()->GetParMesh());
130 subfes1 = src.ParFESpace();
131 parentfes = dst.ParFESpace();
133 *parentfes,
134 src_sm->GetFrom(),
135 src_sm->GetParentElementIDMap(),
136 sub1_to_parent_map_);
137
138 root_gc_ = &parentfes->GroupComm();
139 CommunicateIndicesSet(sub1_to_parent_map_, dst.Size());
140 }
142 {
144 ParSubMesh* dst_sm = static_cast<ParSubMesh*>(dst.ParFESpace()->GetParMesh());
145 subfes1 = dst.ParFESpace();
146 parentfes = src.ParFESpace();
148 *parentfes,
149 dst_sm->GetFrom(),
150 dst_sm->GetParentElementIDMap(),
151 sub1_to_parent_map_);
152 }
153 else
154 {
155 MFEM_ABORT("Trying to do a transfer between GridFunctions but none of them is defined on a SubMesh");
156 }
157}
158
160 ParGridFunction &dst) const
161{
162 if (category_ == TransferCategory::ParentToSubMesh)
163 {
164 // dst = S1^T src
165 src.HostRead();
166 dst.HostWrite(); // dst is fully overwritten
167 for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
168 {
169 real_t s = 1.0;
170 int j = FiniteElementSpace::DecodeDof(sub1_to_parent_map_[i], s);
171 dst(i) = s * src(j);
172 }
173
174 CorrectFaceOrientations(*dst.ParFESpace(), src, dst);
175 }
176 else if (category_ == TransferCategory::SubMeshToParent)
177 {
178 // dst = G S1 src
179 // = G z
180 //
181 // G is identity if the partitioning matches
182
183 src.HostRead();
184 dst.HostReadWrite(); // dst is only partially overwritten
185 for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
186 {
187 real_t s = 1.0;
188 int j = FiniteElementSpace::DecodeDof(sub1_to_parent_map_[i], s);
189 dst(j) = s * src(i);
190 }
191
192 CorrectFaceOrientations(*src.ParFESpace(), src, dst,
193 &sub1_to_parent_map_);
194
195 CommunicateSharedVdofs(dst);
196 }
197 else if (category_ == TransferCategory::SubMeshToSubMesh)
198 {
199 // dst = S2^T G (S1 src (*) S2 dst)
200 //
201 // G is identity if the partitioning matches
202
203 src.HostRead();
204 dst.HostReadWrite();
205
206 z_ = 0.0;
207
208 for (int i = 0; i < sub2_to_parent_map_.Size(); i++)
209 {
210 real_t s = 1.0;
211 int j = FiniteElementSpace::DecodeDof(sub2_to_parent_map_[i], s);
212 z_(j) = s * dst(i);
213 }
214
215 CorrectFaceOrientations(*dst.ParFESpace(), dst, z_,
216 &sub2_to_parent_map_);
217
218 for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
219 {
220 real_t s = 1.0;
221 int j = FiniteElementSpace::DecodeDof(sub1_to_parent_map_[i], s);
222 z_(j) = s * src(i);
223 }
224
225 CorrectFaceOrientations(*src.ParFESpace(), src, z_,
226 &sub1_to_parent_map_);
227
228 CommunicateSharedVdofs(z_);
229
230 for (int i = 0; i < sub2_to_parent_map_.Size(); i++)
231 {
232 real_t s = 1.0;
233 int j = FiniteElementSpace::DecodeDof(sub2_to_parent_map_[i], s);
234 dst(i) = s * z_(j);
235 }
236
237 CorrectFaceOrientations(*dst.ParFESpace(), z_, dst);
238 }
239 else
240 {
241 MFEM_ABORT("unknown TransferCategory: " << category_);
242 }
243}
244
245void ParTransferMap::CommunicateIndicesSet(Array<int> &map, int dst_sz)
246{
247 indices_set_local_.SetSize(dst_sz);
248 indices_set_local_ = 0;
249 for (int i = 0; i < map.Size(); i++)
250 {
251 indices_set_local_[(map[i]>=0)?map[i]:(-map[i]-1)] = 1;
252 }
253 indices_set_global_ = indices_set_local_;
254 root_gc_->Reduce(indices_set_global_, GroupCommunicator::Sum);
255 root_gc_->Bcast(indices_set_global_);
256}
257
258void ParTransferMap::CommunicateSharedVdofs(Vector &f) const
259{
260 // f is usually defined on the root vdofs
261
262 const Table &group_ldof = root_gc_->GroupLDofTable();
263
264 // Identify indices that were only set by other ranks and clear the dof.
265 for (int i = 0; i < group_ldof.Size_of_connections(); i++)
266 {
267 const int j = group_ldof.GetJ()[i];
268 if (indices_set_global_[j] != 0 && indices_set_local_[j] == 0)
269 {
270 f(j) = 0.0;
271 }
272 }
273
274 // TODO: do the reduce only on dofs of interest
275 root_gc_->Reduce<real_t>(f.HostReadWrite(), GroupCommunicator::Sum);
276
277 // Indices that were set from this rank or other ranks have been summed up
278 // and therefore need to be "averaged". Note that this results in the exact
279 // value that is desired.
280 for (int i = 0; i < group_ldof.Size_of_connections(); i++)
281 {
282 const int j = group_ldof.GetJ()[i];
283 if (indices_set_global_[j] != 0)
284 {
285 f(j) /= indices_set_global_[j];
286 }
287 }
288
289 // Indices for dofs that are shared between processors need to be divided by
290 // the whole group size that share this dof.
291 for (int gr = 1; gr < group_ldof.Size(); gr++)
292 {
293 for (int i = 0; i < group_ldof.RowSize(gr); i++)
294 {
295 const int j = group_ldof.GetRow(gr)[i];
296 if (indices_set_global_[j] == 0)
297 {
298 f(j) /= root_gc_->GetGroupTopology().GetGroupSize(gr);
299 }
300 }
301 }
302
303 root_gc_->Bcast<real_t>(f.HostReadWrite());
304}
305
306void
307ParTransferMap::CorrectFaceOrientations(const ParFiniteElementSpace &fes,
308 const Vector &src,
309 Vector &dst,
310 const Array<int> *sub_to_parent_map)
311{
312 const FiniteElementCollection * fec = fes.FEColl();
313
314 ParSubMesh * mesh = dynamic_cast<ParSubMesh*>(fes.GetParMesh());
315
316 const Array<int>& parent_face_ori = mesh->GetParentFaceOrientations();
317
318 if (parent_face_ori.Size() == 0) { return; }
319
320 DofTransformation doftrans(fes.GetVDim(), fes.GetOrdering());
321
322 int dim = mesh->Dimension();
323 bool face = (dim == 3);
324
325 Array<int> vdofs;
326 Array<int> Fo(1);
327 Vector face_vector;
328
329 for (int i = 0; i < (face ? mesh->GetNumFaces() : mesh->GetNE()); i++)
330 {
331 if (parent_face_ori[i] == 0) { continue; }
332
333 Geometry::Type geom = face ? mesh->GetFaceGeometry(i) :
334 mesh->GetElementGeometry(i);
335
336 if (!fec->DofTransformationForGeometry(geom)) { continue; }
337 doftrans.SetDofTransformation(*fec->DofTransformationForGeometry(geom));
338
339 Fo[0] = parent_face_ori[i];
340 doftrans.SetFaceOrientations(Fo);
341
342 if (face)
343 {
344 fes.GetFaceVDofs(i, vdofs);
345 }
346 else
347 {
348 fes.GetElementVDofs(i, vdofs);
349 }
350
351 if (sub_to_parent_map)
352 {
353 src.GetSubVector(vdofs, face_vector);
354 doftrans.TransformPrimal(face_vector);
355 }
356 else
357 {
358 dst.GetSubVector(vdofs, face_vector);
359 doftrans.InvTransformPrimal(face_vector);
360 }
361
362 for (int j = 0; j < vdofs.Size(); j++)
363 {
364 real_t s = 1.0;
365 int k = FiniteElementSpace::DecodeDof(vdofs[j], s);
366
367 if (sub_to_parent_map)
368 {
369 real_t sps = 1.0;
370 int spk = FiniteElementSpace::DecodeDof((*sub_to_parent_map)[k],
371 sps);
372 s *= sps;
373 k = spk;
374 }
375
376 dst[k] = s * face_vector[j];
377 }
378 }
379}
380
381#endif // MFEM_USE_MPI
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
@ GaussLegendre
Open type.
Definition fe_base.hpp:31
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition fe_coll.hpp:225
virtual const StatelessDofTransformation * DofTransformationForGeometry(Geometry::Type GeomType) const
Returns a DoF transformation object compatible with this basis and geometry type.
Definition fe_coll.hpp:67
int GetMapType(int dim) const
Definition fe_coll.cpp:60
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:287
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:725
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition fespace.cpp:316
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Definition fespace.hpp:1017
Table & GroupLDofTable()
Fill-in the returned Table reference to initialize the GroupCommunicator then call Finalize().
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
const GroupTopology & GetGroupTopology()
Get a reference to the associated GroupTopology object.
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
int GetGroupSize(int g) const
Get the number of processors in a group.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
int GetBasisType() const
Definition fe_coll.hpp:373
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6250
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
Definition mesh.cpp:1452
Geometry::Type GetElementGeometry(int i) const
Definition mesh.hpp:1371
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
Abstract parallel finite element space.
Definition pfespace.hpp:29
GroupCommunicator & GroupComm()
Return a reference to the internal GroupCommunicator (on VDofs)
Definition pfespace.hpp:344
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
Class for parallel grid function.
Definition pgridfunc.hpp:33
ParFiniteElementSpace * ParFESpace() const
Class for parallel meshes.
Definition pmesh.hpp:34
Subdomain representation of a topological parent in another ParMesh.
Definition psubmesh.hpp:52
static bool IsParSubMesh(const ParMesh *m)
Check if ParMesh m is a ParSubMesh.
Definition psubmesh.hpp:180
const Array< int > & GetParentElementIDMap() const
Get the parent element id map.
Definition psubmesh.hpp:106
const Array< int > & GetParentFaceOrientations() const
Get the relative face orientations.
Definition psubmesh.hpp:136
SubMesh::From GetFrom() const
Get the From indicator.
Definition psubmesh.hpp:96
void Transfer(const ParGridFunction &src, ParGridFunction &dst) const
Transfer the source ParGridFunction to the destination ParGridFunction.
ParTransferMap(const ParGridFunction &src, const ParGridFunction &dst)
Construct a new ParTransferMap object which transfers degrees of freedom from the source ParGridFunct...
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
int * GetJ()
Definition table.hpp:114
int RowSize(int i) const
Definition table.hpp:108
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition table.cpp:187
int Size() const
Returns the number of TYPE I elements.
Definition table.hpp:92
int Size_of_connections() const
Definition table.hpp:98
Vector data type.
Definition vector.hpp:80
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:478
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:486
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:494
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
int dim
Definition ex24.cpp:53
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.
RT GetRootParent(const T &m)
Identify the root parent of a given SubMesh.
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
RefCoord s[3]