MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
submesh_utils.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_utils.hpp"
13
14namespace mfem
15{
16namespace SubMeshUtils
17{
18int UniqueIndexGenerator::Get(int i, bool &new_index)
19{
20 auto f = idx.find(i);
21 if (f == idx.end())
22 {
23 idx[i] = counter;
24 new_index = true;
25 return counter++;
26 }
27 else
28 {
29 new_index = false;
30 return (*f).second;
31 }
32}
33
34bool ElementHasAttribute(const Element &el, const Array<int> &attributes)
35{
36 for (int a = 0; a < attributes.Size(); a++)
37 {
38 if (el.GetAttribute() == attributes[a])
39 {
40 return true;
41 }
42 }
43 return false;
44}
45
46std::tuple< Array<int>, Array<int> >
47AddElementsToMesh(const Mesh& parent,
48 Mesh& mesh,
49 const Array<int> &attributes,
50 bool from_boundary)
51{
52 Array<int> parent_vertex_ids, parent_element_ids;
53 UniqueIndexGenerator vertex_ids;
54 const int ne = from_boundary ? parent.GetNBE() : parent.GetNE();
55 for (int i = 0; i < ne; i++)
56 {
57 const Element *pel = from_boundary ?
58 parent.GetBdrElement(i) : parent.GetElement(i);
59 if (!ElementHasAttribute(*pel, attributes)) { continue; }
60
61 Array<int> v;
62 pel->GetVertices(v);
63 Array<int> submesh_v(v.Size());
64
65 for (int iv = 0; iv < v.Size(); iv++)
66 {
67 bool new_vertex;
68 int mesh_vertex_id = v[iv];
69 int submesh_vertex_id = vertex_ids.Get(mesh_vertex_id, new_vertex);
70 if (new_vertex)
71 {
72 mesh.AddVertex(parent.GetVertex(mesh_vertex_id));
73 parent_vertex_ids.Append(mesh_vertex_id);
74 }
75 submesh_v[iv] = submesh_vertex_id;
76 }
77
78 Element *el = mesh.NewElement(from_boundary ?
79 parent.GetBdrElementType(i) : parent.GetElementType(i));
80 el->SetVertices(submesh_v);
81 el->SetAttribute(pel->GetAttribute());
82 mesh.AddElement(el);
83 parent_element_ids.Append(i);
84 }
85 return std::tuple<Array<int>, Array<int>>(parent_vertex_ids,
86 parent_element_ids);
87}
88
90 const FiniteElementSpace& parentfes,
91 const SubMesh::From& from,
92 const Array<int>& parent_element_ids,
93 Array<int>& vdof_to_vdof_map)
94{
95 auto *m = subfes.GetMesh();
96 vdof_to_vdof_map.SetSize(subfes.GetVSize());
97
98 const int vdim = parentfes.GetVDim();
99
101 DenseMatrix T;
102 Array<int> z1;
103 for (int i = 0; i < m->GetNE(); i++)
104 {
105 Array<int> parent_vdofs;
106 if (from == SubMesh::From::Domain)
107 {
108 parentfes.GetElementVDofs(parent_element_ids[i], parent_vdofs);
109 }
110 else if (from == SubMesh::From::Boundary)
111 {
112 if (parentfes.IsDGSpace())
113 {
114 MFEM_ASSERT(static_cast<const L2_FECollection*>
115 (parentfes.FEColl())->GetBasisType() == BasisType::GaussLobatto,
116 "Only BasisType::GaussLobatto is supported for L2 spaces");
117
118 auto pm = parentfes.GetMesh();
119
120 const Geometry::Type face_geom =
121 pm->GetBdrElementGeometry(parent_element_ids[i]);
122 int face_info, parent_volel_id;
123 pm->GetBdrElementAdjacentElement(
124 parent_element_ids[i], parent_volel_id, face_info);
125 face_info = Mesh::EncodeFaceInfo(
128 face_geom, Mesh::DecodeFaceInfoOrientation(face_info)));
129 pm->GetLocalFaceTransformation(
130 pm->GetBdrElementType(parent_element_ids[i]),
131 pm->GetElementType(parent_volel_id),
132 Tr.Transf,
133 face_info);
134
135 const FiniteElement *face_el =
136 parentfes.GetTraceElement(parent_element_ids[i], face_geom);
137 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
138 "Nodal Finite Element is required");
139
140 face_el->GetTransferMatrix(*parentfes.GetFE(parent_volel_id),
141 Tr.Transf,
142 T);
143
144 parentfes.GetElementVDofs(parent_volel_id, z1);
145
146 parent_vdofs.SetSize(vdim * T.Height());
147 for (int j = 0; j < T.Height(); j++)
148 {
149 for (int k = 0; k < T.Width(); k++)
150 {
151 if (T(j, k) != 0.0)
152 {
153 for (int vd=0; vd<vdim; vd++)
154 {
155 int sub_vdof = j + T.Height() * vd;
156 int parent_vdof = k + T.Width() * vd;
157 parent_vdofs[sub_vdof] =
158 z1[static_cast<int>(parent_vdof)];
159 }
160 }
161 }
162 }
163 }
164 else
165 {
166 parentfes.GetBdrElementVDofs(parent_element_ids[i], parent_vdofs);
167 }
168 }
169 else
170 {
171 MFEM_ABORT("SubMesh::From type unknown");
172 }
173
174 Array<int> sub_vdofs;
175 subfes.GetElementVDofs(i, sub_vdofs);
176
177 MFEM_ASSERT(parent_vdofs.Size() == sub_vdofs.Size(), "internal error");
178 for (int j = 0; j < parent_vdofs.Size(); j++)
179 {
180 real_t sub_sign = 1.0;
181 int sub_vdof = subfes.DecodeDof(sub_vdofs[j], sub_sign);
182
183 real_t parent_sign = 1.0;
184 int parent_vdof = parentfes.DecodeDof(parent_vdofs[j], parent_sign);
185
186 vdof_to_vdof_map[sub_vdof] =
187 (sub_sign * parent_sign > 0.0) ? parent_vdof : (-1-parent_vdof);
188 }
189 }
190}
191
192Array<int> BuildFaceMap(const Mesh& pm, const Mesh& sm,
193 const Array<int> &parent_element_ids)
194{
195 // TODO: Check if parent is really a parent of mesh
196
197 Array<int> pfids(sm.GetNumFaces());
198 pfids = -1;
199 for (int i = 0; i < sm.GetNE(); i++)
200 {
201 int peid = parent_element_ids[i];
202
203 Array<int> sel_faces, pel_faces, o;
204 if (pm.Dimension() == 2)
205 {
206 sm.GetElementEdges(i, sel_faces, o);
207 pm.GetElementEdges(peid, pel_faces, o);
208 }
209 else
210 {
211 sm.GetElementFaces(i, sel_faces, o);
212 pm.GetElementFaces(peid, pel_faces, o);
213 }
214
215 MFEM_ASSERT(sel_faces.Size() == pel_faces.Size(), "internal error");
216
217 for (int j = 0; j < sel_faces.Size(); j++)
218 {
219 if (pfids[sel_faces[j]] != -1)
220 {
221 MFEM_ASSERT(pfids[sel_faces[j]] == pel_faces[j], "internal error");
222 }
223 pfids[sel_faces[j]] = pel_faces[j];
224 }
225 }
226 return pfids;
227}
228
229} // namespace SubMeshUtils
230} // namespace mfem
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
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Definition element.hpp:58
int GetAttribute() const
Return element's attribute.
Definition element.hpp:55
virtual void SetVertices(const Array< int > &v)=0
Set the indices defining the vertices.
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
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
Definition fespace.cpp:3276
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Definition fespace.hpp:1313
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Definition fespace.hpp:1017
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
Definition fespace.cpp:303
Abstract class for all finite elements.
Definition fe_base.hpp:239
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition fe_base.cpp:119
static int GetInverseOrientation(Type geom_type, int orientation)
Return the inverse of the given orientation for the specified geometry type.
Definition geom.cpp:264
IsoparametricTransformation Transf
Definition eltrans.hpp:467
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
Element * NewElement(int geom)
Definition mesh.cpp:4401
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition mesh.cpp:7316
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
Definition mesh.cpp:7321
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6250
Geometry::Type GetBdrElementGeometry(int i) const
Definition mesh.hpp:1376
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
static int EncodeFaceInfo(int local_face_index, int orientation)
Given local_face_index and orientation, return the corresponding encoded "face info int".
Definition mesh.hpp:1990
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
int AddElement(Element *elem)
Definition mesh.cpp:2021
static int DecodeFaceInfoLocalIndex(int info)
Given a "face info int", return the local face index.
Definition mesh.hpp:1986
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition mesh.cpp:7201
static int DecodeFaceInfoOrientation(int info)
Given a "face info int", return the face orientation.
Definition mesh.hpp:1983
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition mesh.cpp:6985
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
Class for standard nodal finite elements.
Definition fe_base.hpp:715
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
From
Indicator from which part of the parent Mesh the SubMesh is created.
Definition submesh.hpp:47
real_t a
Definition lissajous.cpp:41
std::tuple< Array< int >, Array< int > > AddElementsToMesh(const Mesh &parent, Mesh &mesh, const Array< int > &attributes, bool from_boundary)
Given a Mesh parent and another Mesh mesh using the list of attributes in attributes,...
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.
Array< int > BuildFaceMap(const Mesh &pm, const Mesh &sm, const Array< int > &parent_element_ids)
Given two meshes that have a parent to SubMesh relationship create a face map, using a SubMesh to par...
bool ElementHasAttribute(const Element &el, const Array< int > &attributes)
Given an element el and a list of attributes, determine if that element is in at least one attribute ...
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
Convenience object to create unique indices.
std::unordered_map< int, int > idx
int Get(int i, bool &new_index)
Returns the unique index from an index set.