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 9398 : MFEMCutTransitionSubMesh::validParams()
19 : {
20 9398 : InputParameters params = MFEMSubMesh::validParams();
21 9398 : params += MFEMBlockRestrictable::validParams();
22 18796 : 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 37592 : params.addRequiredParam<BoundaryName>("cut_boundary",
27 : "The boundary associated with the mesh cut.");
28 37592 : params.addRequiredParam<BoundaryName>(
29 : "transition_subdomain_boundary",
30 : "Name to assign boundary of transition subdomain not shared with cut surface.");
31 37592 : 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 28194 : 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 9398 : 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 = getMFEMProblem().getProblemData().myid;
70 23 : int mpi_comm_size = getMFEMProblem().getProblemData().num_procs;
71 :
72 : // First determine face normal based on the first boundary element found on the cut
73 : // to use when determining orientation relative to the cut
74 23 : const mfem::Array<int> & parent_cut_element_id_map = _cut_submesh->GetParentElementIDMap();
75 23 : int rank_with_submesh = -1;
76 23 : if (parent_cut_element_id_map.Size() > 0)
77 : {
78 21 : int reference_face = parent_cut_element_id_map[0];
79 21 : _cut_normal = findFaceNormal(parent_mesh, parent_mesh.GetBdrElementFaceIndex(reference_face));
80 21 : rank_with_submesh = mpi_comm_rank;
81 : }
82 23 : MPI_Allreduce(MPI_IN_PLACE, &rank_with_submesh, 1, MPI_INT, MPI_MAX, getMFEMProblem().getComm());
83 23 : MPI_Bcast(_cut_normal.GetData(),
84 : _cut_normal.Size(),
85 : MFEM_MPI_REAL_T,
86 : rank_with_submesh,
87 : getMFEMProblem().getComm());
88 : // Iterate over all vertices on cut, find elements with those vertices,
89 : // and declare them transition elements if they are on the +ve side of the cut
90 23 : mfem::Array<int> transition_els;
91 23 : std::vector<HYPRE_BigInt> global_cut_vert_ids;
92 23 : mfem::Array<HYPRE_BigInt> gi;
93 23 : parent_mesh.GetGlobalVertexIndices(gi);
94 23 : std::unique_ptr<mfem::Table> vert_to_elem(parent_mesh.GetVertexToElementTable());
95 23 : const mfem::Array<int> & cut_to_parent_vertex_id_map = _cut_submesh->GetParentVertexIDMap();
96 2537 : for (int i = 0; i < _cut_submesh->GetNV(); ++i)
97 : {
98 2514 : int cut_vert = cut_to_parent_vertex_id_map[i];
99 2514 : global_cut_vert_ids.push_back(gi[cut_vert]);
100 2514 : int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
101 2514 : const int * els_adj_to_cut = vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
102 62056 : for (int i = 0; i < ne; i++)
103 : {
104 59542 : const int el_adj_to_cut = els_adj_to_cut[i];
105 108582 : if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
106 49040 : isPositiveSideOfCut(el_adj_to_cut, cut_vert, parent_mesh))
107 24540 : transition_els.Append(el_adj_to_cut);
108 : }
109 : }
110 :
111 : // share cut verts coords or global ids across all procs
112 23 : int n_cut_vertices = global_cut_vert_ids.size();
113 23 : std::vector<int> cut_vert_sizes(mpi_comm_size, 0);
114 23 : MPI_Allgather(
115 : &n_cut_vertices, 1, MPI_INT, &cut_vert_sizes[0], 1, MPI_INT, getMFEMProblem().getComm());
116 : // Make an offset array and total the sizes.
117 23 : std::vector<int> n_vert_offset(mpi_comm_size, 0);
118 35 : for (int i = 1; i < mpi_comm_size; i++)
119 12 : n_vert_offset[i] = n_vert_offset[i - 1] + cut_vert_sizes[i - 1];
120 23 : int global_n_cut_vertices = 0;
121 58 : for (int i = 0; i < mpi_comm_size; i++)
122 35 : global_n_cut_vertices += cut_vert_sizes[i];
123 :
124 : // Gather the queries to all ranks.
125 23 : std::vector<HYPRE_BigInt> all_cut_verts(global_n_cut_vertices, 0);
126 23 : MPI_Allgatherv(&global_cut_vert_ids[0],
127 : n_cut_vertices,
128 : HYPRE_MPI_BIG_INT,
129 : &all_cut_verts[0],
130 : &cut_vert_sizes[0],
131 : &n_vert_offset[0],
132 : HYPRE_MPI_BIG_INT,
133 : getMFEMProblem().getComm());
134 :
135 : // Detect shared vertices and add corresponding elements
136 35 : for (int g = 1, sv = 0; g < parent_mesh.GetNGroups(); g++)
137 : {
138 6392 : for (int gv = 0; gv < parent_mesh.GroupNVertices(g); gv++, sv++)
139 : {
140 : // all elements touching this shared vertex should be updated
141 6380 : int cut_vert = parent_mesh.GroupVertex(g, gv);
142 941008 : for (std::size_t i = 0; i < all_cut_verts.size(); i += 1)
143 : {
144 934628 : if (gi[cut_vert] == all_cut_verts[i]) // check if shared vertex is on the cut plane
145 : {
146 208 : int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
147 : const int * els_adj_to_cut =
148 208 : vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
149 2880 : for (int i = 0; i < ne; i++)
150 : {
151 2672 : const int el_adj_to_cut = els_adj_to_cut[i];
152 5104 : if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
153 2432 : isPositiveSideOfCut(el_adj_to_cut, cut_vert, parent_mesh))
154 1208 : transition_els.Append(el_adj_to_cut);
155 : }
156 : }
157 : }
158 : }
159 : }
160 :
161 23 : transition_els.Sort();
162 23 : transition_els.Unique();
163 :
164 23 : setAttributes(parent_mesh, transition_els);
165 23 : }
166 :
167 : mfem::Vector
168 21 : MFEMCutTransitionSubMesh::findFaceNormal(const mfem::ParMesh & mesh, const int & face)
169 : {
170 : mooseAssert(mesh.SpaceDimension() == 3,
171 : "MFEMCutTransitionSubMesh only works in 3-dimensional meshes!");
172 21 : mfem::Vector normal;
173 21 : mfem::Array<int> face_verts;
174 21 : std::vector<mfem::Vector> v;
175 21 : mesh.GetFaceVertices(face, face_verts);
176 :
177 : // First we get the coordinates of 3 vertices on the face
178 84 : for (auto vtx : face_verts)
179 : {
180 63 : mfem::Vector vtx_coords(3);
181 252 : for (int j = 0; j < 3; ++j)
182 189 : vtx_coords[j] = mesh.GetVertex(vtx)[j];
183 63 : v.push_back(vtx_coords);
184 63 : }
185 :
186 : // Now we find the unit vector normal to the face
187 21 : v[0] -= v[1];
188 21 : v[1] -= v[2];
189 21 : v[0].cross3D(v[1], normal);
190 21 : normal /= normal.Norml2();
191 42 : return normal;
192 21 : }
193 :
194 : bool
195 51472 : MFEMCutTransitionSubMesh::isPositiveSideOfCut(const int & el,
196 : const int & el_vertex_on_cut,
197 : mfem::ParMesh & parent_mesh)
198 : {
199 51472 : const int sdim = parent_mesh.SpaceDimension();
200 51472 : mfem::Vector el_center(3);
201 51472 : parent_mesh.GetElementCenter(el, el_center);
202 51472 : mfem::Vector vertex_coords(parent_mesh.GetVertex(el_vertex_on_cut), sdim);
203 51472 : mfem::Vector relative_center(sdim);
204 205888 : for (int j = 0; j < sdim; j++)
205 154416 : relative_center[j] = el_center[j] - vertex_coords[j];
206 102944 : return _cut_normal * relative_center > 0;
207 51472 : }
208 :
209 : void
210 23 : MFEMCutTransitionSubMesh::setAttributes(mfem::ParMesh & parent_mesh,
211 : mfem::Array<int> & transition_els)
212 : {
213 : // Generate a set of local new attributes for transition region elements
214 23 : const int old_max_attr = parent_mesh.attributes.Max();
215 23 : mfem::Array<int> new_attrs(old_max_attr);
216 23 : new_attrs = -1;
217 12337 : for (auto const & transition_el : transition_els)
218 : {
219 12314 : const int old_attr = parent_mesh.GetAttribute(transition_el);
220 12314 : new_attrs[old_attr - 1] = old_max_attr + old_attr;
221 : // Set the new attribute IDs for transition region elements
222 12314 : parent_mesh.SetAttribute(transition_el, new_attrs[old_attr - 1]);
223 : }
224 :
225 : // Distribute local attribute IDs to other MPI ranks to create set of new
226 : // global attribute IDs for the transition region.
227 23 : MPI_Allreduce(
228 : MPI_IN_PLACE, new_attrs, old_max_attr, MPI_INT, MPI_MAX, getMFEMProblem().getComm());
229 :
230 23 : mfem::AttributeSets & attr_sets = parent_mesh.attribute_sets;
231 23 : mfem::AttributeSets & bdr_attr_sets = parent_mesh.bdr_attribute_sets;
232 : // Create attribute set labelling the newly created transition region on one side of the cut
233 23 : attr_sets.CreateAttributeSet(_transition_subdomain);
234 : // Create attribute set labelling the entire closed geometry
235 23 : attr_sets.SetAttributeSet(_closed_subdomain, getSubdomainAttributes());
236 : // Add the new domain attributes to new attribute sets
237 23 : const std::set<std::string> attr_set_names = attr_sets.GetAttributeSetNames();
238 92 : for (int old_attr = 1; old_attr <= old_max_attr; ++old_attr)
239 : {
240 69 : int new_attr = new_attrs[old_attr - 1];
241 : // Add attributes only if they're transition region attributes
242 69 : if (new_attr != -1)
243 : {
244 31 : attr_sets.AddToAttributeSet(_transition_subdomain, new_attr);
245 186 : for (auto const & attr_set_name : attr_set_names)
246 : {
247 155 : const mfem::Array<int> & attr_set = attr_sets.GetAttributeSet(attr_set_name);
248 : // If attribute set has the old attribute of the transition region, add the new one
249 155 : if (attr_set.Find(old_attr) != -1)
250 62 : attr_sets.AddToAttributeSet(attr_set_name, new_attr);
251 : }
252 : }
253 : }
254 :
255 : // Create boundary attribute set labelling the exterior of the newly created
256 : // transition region, excluding the cut
257 23 : bdr_attr_sets.SetAttributeSet(_transition_subdomain_boundary,
258 46 : mfem::Array<int>({parent_mesh.bdr_attributes.Max() + 1}));
259 :
260 23 : parent_mesh.SetAttributes();
261 23 : }
262 :
263 : bool
264 62214 : MFEMCutTransitionSubMesh::isInDomain(const int & element,
265 : const mfem::Array<int> & subdomains,
266 : const mfem::ParMesh & mesh)
267 : {
268 : // element<0 for ghost elements
269 62214 : if (element < 0)
270 0 : return true;
271 :
272 78590 : for (const auto & subdomain : subdomains)
273 67848 : if (mesh.GetAttribute(element) == subdomain)
274 51472 : return true;
275 10742 : return false;
276 : }
277 :
278 : #endif
|
