MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
transfermap.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
12#include "submesh.hpp"
13#include "transfermap.hpp"
14#include "submesh_utils.hpp"
15
16using namespace mfem;
17
19 const GridFunction &dst)
20{
21 const FiniteElementSpace *parentfes = nullptr, *subfes1 = nullptr,
22 *subfes2 = nullptr;
23
24 if (SubMesh::IsSubMesh(src.FESpace()->GetMesh()) &&
26 {
27 SubMesh* src_sm = static_cast<SubMesh*>(src.FESpace()->GetMesh());
28 SubMesh* dst_sm = static_cast<SubMesh*>(dst.FESpace()->GetMesh());
29
30 // There is no immediate relation and both src and dst come from a
31 // SubMesh, check if they have an equivalent root parent.
32 if (SubMeshUtils::GetRootParent(*src_sm) !=
34 {
35 MFEM_ABORT("Can't find a relation between the two GridFunctions");
36 }
37
39
40 {
41 Mesh * parent_mesh =
42 const_cast<Mesh *>(SubMeshUtils::GetRootParent(*src_sm));
43
44 int parent_dim = parent_mesh->Dimension();
45 int src_sm_dim = src_sm->Dimension();
46 int dst_sm_dim = dst_sm->Dimension();
47
48 bool root_fes_reset = false;
49 if (src_sm_dim == parent_dim - 1 && dst_sm_dim == parent_dim - 1)
50 {
51 const FiniteElementSpace *src_fes = src.FESpace();
52 const FiniteElementSpace *dst_fes = dst.FESpace();
53
54 const FiniteElementCollection *src_fec = src_fes->FEColl();
55 const FiniteElementCollection *dst_fec = dst_fes->FEColl();
56
57 const L2_FECollection *src_l2_fec =
58 dynamic_cast<const L2_FECollection*>(src_fec);
59 const L2_FECollection *dst_l2_fec =
60 dynamic_cast<const L2_FECollection*>(dst_fec);
61
62 if (src_l2_fec != NULL && dst_l2_fec != NULL)
63 {
64 // Source and destination are both lower dimension L2 spaces.
65 // Transfer them as the trace of an RT space if possible.
66
67 int src_mt = src_fec->GetMapType(src_sm_dim);
68 int dst_mt = dst_fec->GetMapType(dst_sm_dim);
69
70 int src_bt = src_l2_fec->GetBasisType();
71 int dst_bt = dst_l2_fec->GetBasisType();
72
73 int src_p = src_fec->GetOrder();
74 int dst_p = dst_fec->GetOrder();
75
76 if (src_mt == FiniteElement::INTEGRAL &&
77 dst_mt == FiniteElement::INTEGRAL &&
78 src_bt == BasisType::GaussLegendre &&
79 dst_bt == BasisType::GaussLegendre &&
80 src_p == dst_p)
81 {
82 // The subspaces are consistent with the trace of an RT space
83 root_fec_.reset(new RT_FECollection(src_p, parent_dim));
84 root_fes_.reset(new FiniteElementSpace(
85 const_cast<Mesh *>(
87 root_fec_.get()));
88 root_fes_reset = true;
89 }
90 }
91 }
92
93 if (!root_fes_reset)
94 {
95 root_fes_.reset(new FiniteElementSpace(
96 *src.FESpace(),
97 const_cast<Mesh *>(
99 }
100 }
101
102 subfes1 = src.FESpace();
103 subfes2 = dst.FESpace();
104
106 *root_fes_,
107 src_sm->GetFrom(),
108 src_sm->GetParentElementIDMap(),
109 sub1_to_parent_map_);
110
112 *root_fes_,
113 dst_sm->GetFrom(),
114 dst_sm->GetParentElementIDMap(),
115 sub2_to_parent_map_);
116
117 z_.SetSize(root_fes_->GetVSize());
118 }
119 else if (SubMesh::IsSubMesh(src.FESpace()->GetMesh()))
120 {
122 SubMesh* src_sm = static_cast<SubMesh*>(src.FESpace()->GetMesh());
123 subfes1 = src.FESpace();
124 parentfes = dst.FESpace();
126 *parentfes,
127 src_sm->GetFrom(),
128 src_sm->GetParentElementIDMap(),
129 sub1_to_parent_map_);
130 }
131 else if (SubMesh::IsSubMesh(dst.FESpace()->GetMesh()))
132 {
134 SubMesh* dst_sm = static_cast<SubMesh*>(dst.FESpace()->GetMesh());
135 subfes1 = dst.FESpace();
136 parentfes = src.FESpace();
138 *parentfes,
139 dst_sm->GetFrom(),
140 dst_sm->GetParentElementIDMap(),
141 sub1_to_parent_map_);
142 }
143 else
144 {
145 MFEM_ABORT("Trying to do a transfer between GridFunctions but none of them is defined on a SubMesh");
146 }
147}
148
150 GridFunction &dst) const
151{
152 if (category_ == TransferCategory::ParentToSubMesh)
153 {
154 // dst = S1^T src
155 src.HostRead();
156 dst.HostWrite(); // dst is fully overwritten
157 for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
158 {
159 real_t s = 1.0;
160 int j = FiniteElementSpace::DecodeDof(sub1_to_parent_map_[i], s);
161 dst(i) = s * src(j);
162 }
163
164 CorrectFaceOrientations(*dst.FESpace(), src, dst);
165 }
166 else if (category_ == TransferCategory::SubMeshToParent)
167 {
168 // dst = G S1 src
169 // = G z
170 //
171 // G is identity if the partitioning matches
172
173 src.HostRead();
174 dst.HostReadWrite(); // dst is only partially overwritten
175 for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
176 {
177 real_t s = 1.0;
178 int j = FiniteElementSpace::DecodeDof(sub1_to_parent_map_[i], s);
179 dst(j) = s * src(i);
180 }
181
182 CorrectFaceOrientations(*src.FESpace(), src, dst,
183 &sub1_to_parent_map_);
184 }
185 else if (category_ == TransferCategory::SubMeshToSubMesh)
186 {
187 // dst = S2^T G (S1 src (*) S2 dst)
188 //
189 // G is identity if the partitioning matches
190
191 src.HostRead();
192 dst.HostReadWrite();
193
194 z_ = 0.0;
195
196 for (int i = 0; i < sub2_to_parent_map_.Size(); i++)
197 {
198 real_t s = 1.0;
199 int j = FiniteElementSpace::DecodeDof(sub2_to_parent_map_[i], s);
200 z_(j) = s * dst(i);
201 }
202
203 CorrectFaceOrientations(*dst.FESpace(), dst, z_,
204 &sub2_to_parent_map_);
205
206 for (int i = 0; i < sub1_to_parent_map_.Size(); i++)
207 {
208 real_t s = 1.0;
209 int j = FiniteElementSpace::DecodeDof(sub1_to_parent_map_[i], s);
210 z_(j) = s * src(i);
211 }
212
213 CorrectFaceOrientations(*src.FESpace(), src, z_,
214 &sub1_to_parent_map_);
215
216 for (int i = 0; i < sub2_to_parent_map_.Size(); i++)
217 {
218 real_t s = 1.0;
219 int j = FiniteElementSpace::DecodeDof(sub2_to_parent_map_[i], s);
220 dst(i) = s * z_(j);
221 }
222
223 CorrectFaceOrientations(*dst.FESpace(), z_, dst);
224 }
225 else
226 {
227 MFEM_ABORT("unknown TransferCategory: " << category_);
228 }
229}
230
231void TransferMap::CorrectFaceOrientations(const FiniteElementSpace &fes,
232 const Vector &src,
233 Vector &dst,
234 const Array<int> *sub_to_parent_map)
235{
236 const FiniteElementCollection * fec = fes.FEColl();
237
238 SubMesh * mesh = dynamic_cast<SubMesh*>(fes.GetMesh());
239
240 const Array<int>& parent_face_ori = mesh->GetParentFaceOrientations();
241
242 if (parent_face_ori.Size() == 0) { return; }
243
244 DofTransformation doftrans(fes.GetVDim(), fes.GetOrdering());
245
246 int dim = mesh->Dimension();
247 bool face = (dim == 3);
248
249 Array<int> vdofs;
250 Array<int> Fo(1);
251 Vector face_vector;
252
253 for (int i = 0; i < (face ? mesh->GetNumFaces() : mesh->GetNE()); i++)
254 {
255 if (parent_face_ori[i] == 0) { continue; }
256
257 Geometry::Type geom = face ? mesh->GetFaceGeometry(i) :
258 mesh->GetElementGeometry(i);
259
260 if (!fec->DofTransformationForGeometry(geom)) { continue; }
261 doftrans.SetDofTransformation(*fec->DofTransformationForGeometry(geom));
262
263 Fo[0] = parent_face_ori[i];
264 doftrans.SetFaceOrientations(Fo);
265
266 if (face)
267 {
268 fes.GetFaceVDofs(i, vdofs);
269 }
270 else
271 {
272 fes.GetElementVDofs(i, vdofs);
273 }
274
275 if (sub_to_parent_map)
276 {
277 src.GetSubVector(vdofs, face_vector);
278 doftrans.TransformPrimal(face_vector);
279 }
280 else
281 {
282 dst.GetSubVector(vdofs, face_vector);
283 doftrans.InvTransformPrimal(face_vector);
284 }
285
286 for (int j = 0; j < vdofs.Size(); j++)
287 {
288 real_t s = 1.0;
289 int k = FiniteElementSpace::DecodeDof(vdofs[j], s);
290
291 if (sub_to_parent_map)
292 {
293 real_t sps = 1.0;
294 int spk = FiniteElementSpace::DecodeDof((*sub_to_parent_map)[k],
295 sps);
296 s *= sps;
297 k = spk;
298 }
299
300 dst[k] = s * face_vector[j];
301 }
302 }
303}
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
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
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
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
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
int GetBasisType() const
Definition fe_coll.hpp:373
Mesh data type.
Definition mesh.hpp:56
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
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Subdomain representation of a topological parent in another Mesh.
Definition submesh.hpp:43
const Array< int > & GetParentElementIDMap() const
Get the parent element id map.
Definition submesh.hpp:108
From GetFrom() const
Get the From indicator.
Definition submesh.hpp:98
const Array< int > & GetParentFaceOrientations() const
Get the relative face orientations.
Definition submesh.hpp:128
static bool IsSubMesh(const Mesh *m)
Check if Mesh m is a SubMesh.
Definition submesh.hpp:172
void Transfer(const GridFunction &src, GridFunction &dst) const
Transfer the source GridFunction to the destination GridFunction.
TransferMap(const GridFunction &src, const GridFunction &dst)
Construct a new TransferMap object which transfers degrees of freedom from the source GridFunction to...
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
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
RefCoord s[3]