Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #ifdef MOOSE_MFEM_ENABLED
11 :
12 : #include "MFEMCutTransitionSubMesh.h"
13 : #include "MFEMProblem.h"
14 :
15 : registerMooseObject("MooseApp", MFEMCutTransitionSubMesh);
16 :
17 : InputParameters
18 8676 : MFEMCutTransitionSubMesh::validParams()
19 : {
20 8676 : InputParameters params = MFEMSubMesh::validParams();
21 8676 : params += MFEMBlockRestrictable::validParams();
22 17352 : params.addClassDescription(
23 : "Class to construct an MFEMSubMesh formed from the set of elements that have least one "
24 : "vertex on the specified cut plane, that lie on one side of the plane, "
25 : "and that are restricted to the set of user-specified subdomains.");
26 34704 : params.addRequiredParam<BoundaryName>("cut_boundary",
27 : "The boundary associated with the mesh cut.");
28 34704 : params.addRequiredParam<BoundaryName>(
29 : "transition_subdomain_boundary",
30 : "Name to assign boundary of transition subdomain not shared with cut surface.");
31 34704 : params.addRequiredParam<SubdomainName>(
32 : "transition_subdomain",
33 : "The name of the subdomain to be created on the mesh comprised of the set of elements "
34 : "adjacent to the cut surface on one side.");
35 26028 : params.addRequiredParam<SubdomainName>(
36 : "closed_subdomain",
37 : "The name of the subdomain attribute to be created comprised of the set of all elements "
38 : "of the closed geometry, including the new transition region.");
39 8676 : return params;
40 0 : }
41 :
42 23 : MFEMCutTransitionSubMesh::MFEMCutTransitionSubMesh(const InputParameters & parameters)
43 : : MFEMSubMesh(parameters),
44 23 : MFEMBlockRestrictable(parameters, getMFEMProblem().mesh().getMFEMParMesh()),
45 23 : _cut_boundary(getParam<BoundaryName>("cut_boundary")),
46 23 : _cut_submesh(std::make_shared<mfem::ParSubMesh>(mfem::ParSubMesh::CreateFromBoundary(
47 23 : getMFEMProblem().mesh().getMFEMParMesh(),
48 23 : getMesh().bdr_attribute_sets.GetAttributeSet(_cut_boundary)))),
49 46 : _transition_subdomain_boundary(getParam<BoundaryName>("transition_subdomain_boundary")),
50 46 : _transition_subdomain(getParam<SubdomainName>("transition_subdomain")),
51 46 : _closed_subdomain(getParam<SubdomainName>("closed_subdomain")),
52 69 : _cut_normal(3)
53 : {
54 23 : }
55 :
56 : void
57 23 : MFEMCutTransitionSubMesh::buildSubMesh()
58 : {
59 23 : labelMesh(getMFEMProblem().mesh().getMFEMParMesh());
60 46 : _submesh = std::make_shared<mfem::ParSubMesh>(mfem::ParSubMesh::CreateFromDomain(
61 46 : getMesh(), getMesh().attribute_sets.GetAttributeSet(_transition_subdomain)));
62 23 : _submesh->attribute_sets.attr_sets = getMesh().attribute_sets.attr_sets;
63 23 : _submesh->bdr_attribute_sets.attr_sets = getMesh().bdr_attribute_sets.attr_sets;
64 23 : }
65 :
66 : void
67 23 : MFEMCutTransitionSubMesh::labelMesh(mfem::ParMesh & parent_mesh)
68 : {
69 23 : int mpi_comm_rank = 0;
70 23 : int mpi_comm_size = 1;
71 23 : MPI_Comm_rank(getMFEMProblem().mesh().getMFEMParMesh().GetComm(), &mpi_comm_rank);
72 23 : MPI_Comm_size(getMFEMProblem().mesh().getMFEMParMesh().GetComm(), &mpi_comm_size);
73 :
74 : // First determine face normal based on the first boundary element found on the cut
75 : // to use when determining orientation relative to the cut
76 23 : const mfem::Array<int> & parent_cut_element_id_map = _cut_submesh->GetParentElementIDMap();
77 23 : int rank_with_submesh = -1;
78 23 : if (parent_cut_element_id_map.Size() > 0)
79 : {
80 21 : int reference_face = parent_cut_element_id_map[0];
81 21 : _cut_normal = findFaceNormal(parent_mesh, parent_mesh.GetBdrElementFaceIndex(reference_face));
82 21 : rank_with_submesh = mpi_comm_rank;
83 : }
84 23 : MPI_Allreduce(MPI_IN_PLACE,
85 : &rank_with_submesh,
86 : 1,
87 : MPI_INT,
88 : MPI_MAX,
89 : getMFEMProblem().mesh().getMFEMParMesh().GetComm());
90 23 : MPI_Bcast(_cut_normal.GetData(),
91 : _cut_normal.Size(),
92 : MPI_DOUBLE,
93 : rank_with_submesh,
94 : getMFEMProblem().mesh().getMFEMParMesh().GetComm());
95 : // Iterate over all vertices on cut, find elements with those vertices,
96 : // and declare them transition elements if they are on the +ve side of the cut
97 23 : mfem::Array<int> transition_els;
98 23 : std::vector<HYPRE_BigInt> global_cut_vert_ids;
99 23 : mfem::Array<HYPRE_BigInt> gi;
100 23 : parent_mesh.GetGlobalVertexIndices(gi);
101 23 : std::unique_ptr<mfem::Table> vert_to_elem(parent_mesh.GetVertexToElementTable());
102 23 : const mfem::Array<int> & cut_to_parent_vertex_id_map = _cut_submesh->GetParentVertexIDMap();
103 2537 : for (int i = 0; i < _cut_submesh->GetNV(); ++i)
104 : {
105 2514 : int cut_vert = cut_to_parent_vertex_id_map[i];
106 2514 : global_cut_vert_ids.push_back(gi[cut_vert]);
107 2514 : int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
108 2514 : const int * els_adj_to_cut = vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
109 62056 : for (int i = 0; i < ne; i++)
110 : {
111 59542 : const int el_adj_to_cut = els_adj_to_cut[i];
112 108582 : if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
113 49040 : sideOfCut(el_adj_to_cut, cut_vert, parent_mesh) == 1)
114 24540 : transition_els.Append(el_adj_to_cut);
115 : }
116 : }
117 :
118 : // share cut verts coords or global ids across all procs
119 23 : int n_cut_vertices = global_cut_vert_ids.size();
120 23 : std::vector<int> cut_vert_sizes(mpi_comm_size, 0);
121 23 : MPI_Allgather(&n_cut_vertices,
122 : 1,
123 : MPI_INT,
124 : &cut_vert_sizes[0],
125 : 1,
126 : MPI_INT,
127 : getMFEMProblem().mesh().getMFEMParMesh().GetComm());
128 : // Make an offset array and total the sizes.
129 23 : std::vector<int> n_vert_offset(mpi_comm_size, 0);
130 35 : for (int i = 1; i < mpi_comm_size; i++)
131 12 : n_vert_offset[i] = n_vert_offset[i - 1] + cut_vert_sizes[i - 1];
132 23 : int global_n_cut_vertices = 0;
133 58 : for (int i = 0; i < mpi_comm_size; i++)
134 35 : global_n_cut_vertices += cut_vert_sizes[i];
135 :
136 : // Gather the queries to all ranks.
137 23 : std::vector<HYPRE_BigInt> all_cut_verts(global_n_cut_vertices, 0);
138 23 : MPI_Allgatherv(&global_cut_vert_ids[0],
139 : n_cut_vertices,
140 : HYPRE_MPI_BIG_INT,
141 : &all_cut_verts[0],
142 : &cut_vert_sizes[0],
143 : &n_vert_offset[0],
144 : HYPRE_MPI_BIG_INT,
145 : getMFEMProblem().mesh().getMFEMParMesh().GetComm());
146 :
147 : // Detect shared vertices and add corresponding elements
148 35 : for (int g = 1, sv = 0; g < parent_mesh.GetNGroups(); g++)
149 : {
150 6392 : for (int gv = 0; gv < parent_mesh.GroupNVertices(g); gv++, sv++)
151 : {
152 : // all elements touching this shared vertex should be updated
153 6380 : int cut_vert = parent_mesh.GroupVertex(g, gv);
154 941008 : for (std::size_t i = 0; i < all_cut_verts.size(); i += 1)
155 : {
156 934628 : if (gi[cut_vert] == all_cut_verts[i]) // check if shared vertex is on the cut plane
157 : {
158 208 : int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
159 : const int * els_adj_to_cut =
160 208 : vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
161 2880 : for (int i = 0; i < ne; i++)
162 : {
163 2672 : const int el_adj_to_cut = els_adj_to_cut[i];
164 5104 : if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
165 2432 : sideOfCut(el_adj_to_cut, cut_vert, parent_mesh) == 1)
166 1208 : transition_els.Append(el_adj_to_cut);
167 : }
168 : }
169 : }
170 : }
171 : }
172 :
173 23 : transition_els.Sort();
174 23 : transition_els.Unique();
175 :
176 23 : setAttributes(parent_mesh, transition_els);
177 23 : }
178 :
179 : mfem::Vector
180 21 : MFEMCutTransitionSubMesh::findFaceNormal(const mfem::ParMesh & mesh, const int & face)
181 : {
182 : MFEM_ASSERT(mesh->SpatialDimension() == 3,
183 : "MFEMCutTransitionSubMesh only works in 3-dimensional meshes!");
184 21 : mfem::Vector normal;
185 21 : mfem::Array<int> face_verts;
186 21 : std::vector<mfem::Vector> v;
187 21 : mesh.GetFaceVertices(face, face_verts);
188 :
189 : // First we get the coordinates of 3 vertices on the face
190 84 : for (auto vtx : face_verts)
191 : {
192 63 : mfem::Vector vtx_coords(3);
193 252 : for (int j = 0; j < 3; ++j)
194 189 : vtx_coords[j] = mesh.GetVertex(vtx)[j];
195 63 : v.push_back(vtx_coords);
196 63 : }
197 :
198 : // Now we find the unit vector normal to the face
199 21 : v[0] -= v[1];
200 21 : v[1] -= v[2];
201 21 : v[0].cross3D(v[1], normal);
202 21 : normal /= normal.Norml2();
203 42 : return normal;
204 21 : }
205 :
206 : int
207 51472 : MFEMCutTransitionSubMesh::sideOfCut(const int & el,
208 : const int & el_vertex_on_cut,
209 : mfem::ParMesh & parent_mesh)
210 : {
211 51472 : const int sdim = parent_mesh.SpaceDimension();
212 51472 : mfem::Vector el_center(3);
213 51472 : parent_mesh.GetElementCenter(el, el_center);
214 51472 : mfem::Vector vertex_coords(parent_mesh.GetVertex(el_vertex_on_cut), sdim);
215 51472 : mfem::Vector relative_center(sdim);
216 205888 : for (int j = 0; j < sdim; j++)
217 154416 : relative_center[j] = el_center[j] - vertex_coords[j];
218 51472 : double side = _cut_normal * relative_center;
219 51472 : if (side > 0)
220 25748 : return 1;
221 : else
222 25724 : return -1;
223 51472 : }
224 :
225 : void
226 23 : MFEMCutTransitionSubMesh::setAttributes(mfem::ParMesh & parent_mesh,
227 : mfem::Array<int> & transition_els)
228 : {
229 : // Generate a set of local new attributes for transition region elements
230 23 : const int old_max_attr = parent_mesh.attributes.Max();
231 23 : mfem::Array<int> new_attrs(old_max_attr);
232 92 : for (int i = 0; i < new_attrs.Size(); ++i)
233 69 : new_attrs[i] = i + 1;
234 12337 : for (int i = 0; i < transition_els.Size(); ++i)
235 : {
236 12314 : const auto & transition_el = transition_els[i];
237 12314 : const int old_attr = parent_mesh.GetAttribute(transition_el);
238 12314 : new_attrs[old_attr - 1] = old_max_attr + old_attr % old_max_attr + 1;
239 : }
240 :
241 : // Distribute local attribute IDs to other MPI ranks to create set of new
242 : // global attribute IDs for the transition region.
243 23 : mfem::Array<int> global_new_attrs(old_max_attr);
244 23 : global_new_attrs = -1;
245 23 : MPI_Allreduce(new_attrs,
246 : global_new_attrs,
247 : old_max_attr,
248 : MPI_INT,
249 : MPI_MAX,
250 : getMFEMProblem().mesh().getMFEMParMesh().GetComm());
251 :
252 : // Set the new attribute IDs for transition region elements
253 12337 : for (int i = 0; i < transition_els.Size(); ++i)
254 : {
255 12314 : const auto & transition_el = transition_els[i];
256 12314 : const int old_attr = parent_mesh.GetAttribute(transition_el);
257 12314 : if (old_attr < old_max_attr + 1)
258 12314 : parent_mesh.SetAttribute(transition_el, global_new_attrs[old_attr - 1]);
259 : }
260 :
261 23 : mfem::AttributeSets & attr_sets = parent_mesh.attribute_sets;
262 23 : mfem::AttributeSets & bdr_attr_sets = parent_mesh.bdr_attribute_sets;
263 : // Create attribute set labelling the newly created transition region
264 : // on one side of the cut
265 23 : attr_sets.CreateAttributeSet(_transition_subdomain);
266 : // Create attribute set labelling the entire closed geometry
267 23 : attr_sets.SetAttributeSet(_closed_subdomain, getSubdomainAttributes());
268 : // Add the new domain attributes to new attribute sets
269 23 : const std::set<std::string> attr_set_names = attr_sets.GetAttributeSetNames();
270 92 : for (int old_attr = 1; old_attr < global_new_attrs.Size() + 1; ++old_attr)
271 : {
272 69 : int new_attr = global_new_attrs[old_attr - 1];
273 : // Add attributes only if they're transition region attributes
274 69 : if (new_attr > old_max_attr)
275 : {
276 31 : attr_sets.AddToAttributeSet(_transition_subdomain, new_attr);
277 186 : for (auto const & attr_set_name : attr_set_names)
278 : {
279 155 : const mfem::Array<int> & attr_set = attr_sets.GetAttributeSet(attr_set_name);
280 : // If attribute set has the old attribute of the transition region, add the new one
281 155 : if (attr_set.Find(old_attr) != -1)
282 62 : attr_sets.AddToAttributeSet(attr_set_name, new_attr);
283 : }
284 : }
285 : }
286 :
287 : // Create boundary attribute set labelling the exterior of the newly created
288 : // transition region, excluding the cut
289 23 : bdr_attr_sets.CreateAttributeSet(_transition_subdomain_boundary);
290 23 : bdr_attr_sets.AddToAttributeSet(_transition_subdomain_boundary,
291 23 : parent_mesh.bdr_attributes.Max() + 1);
292 :
293 23 : parent_mesh.SetAttributes();
294 23 : }
295 :
296 : bool
297 62214 : MFEMCutTransitionSubMesh::isInDomain(const int & element,
298 : const mfem::Array<int> & subdomains,
299 : const mfem::ParMesh & mesh)
300 : {
301 : // element<0 for ghost elements
302 62214 : if (element < 0)
303 0 : return true;
304 :
305 78590 : for (const auto & subdomain : subdomains)
306 67848 : if (mesh.GetAttribute(element) == subdomain)
307 51472 : return true;
308 10742 : return false;
309 : }
310 :
311 : #endif
|
