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 2224 : MFEMCutTransitionSubMesh::validParams()
19 : {
20 2224 : InputParameters params = MFEMSubMesh::validParams();
21 2224 : params += MFEMBlockRestrictable::validParams();
22 4448 : 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 8896 : params.addRequiredParam<BoundaryName>("cut_boundary",
27 : "The boundary associated with the mesh cut.");
28 8896 : params.addRequiredParam<BoundaryName>(
29 : "transition_subdomain_boundary",
30 : "Name to assign boundary of transition subdomain not shared with cut surface.");
31 8896 : 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 6672 : 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 2224 : return params;
40 0 : }
41 :
42 63 : MFEMCutTransitionSubMesh::MFEMCutTransitionSubMesh(const InputParameters & parameters)
43 : : MFEMSubMesh(parameters),
44 63 : MFEMBlockRestrictable(parameters, getMFEMProblem().mesh().getMFEMParMesh()),
45 63 : _cut_boundary(getParam<BoundaryName>("cut_boundary")),
46 63 : _cut_submesh(std::make_shared<mfem::ParSubMesh>(mfem::ParSubMesh::CreateFromBoundary(
47 63 : getMesh(), getMesh().bdr_attribute_sets.GetAttributeSet(_cut_boundary)))),
48 126 : _transition_subdomain_boundary(getParam<BoundaryName>("transition_subdomain_boundary")),
49 126 : _transition_subdomain(getParam<SubdomainName>("transition_subdomain")),
50 126 : _closed_subdomain(getParam<SubdomainName>("closed_subdomain")),
51 189 : _cut_normal(3)
52 : {
53 63 : }
54 :
55 : void
56 63 : MFEMCutTransitionSubMesh::buildSubMesh()
57 : {
58 63 : labelMesh(const_cast<mfem::ParMesh &>(getMesh()));
59 126 : _submesh = std::make_shared<mfem::ParSubMesh>(mfem::ParSubMesh::CreateFromDomain(
60 126 : getMesh(), getMesh().attribute_sets.GetAttributeSet(_transition_subdomain)));
61 63 : _submesh->attribute_sets.attr_sets = getMesh().attribute_sets.attr_sets;
62 63 : _submesh->bdr_attribute_sets.attr_sets = getMesh().bdr_attribute_sets.attr_sets;
63 : // Create boundary attribute set labelling the exterior of the newly created
64 : // transition region, excluding the cut
65 63 : _submesh->bdr_attribute_sets.SetAttributeSet(
66 63 : _transition_subdomain_boundary, mfem::Array<int>({getMesh().bdr_attributes.Max() + 1}));
67 63 : }
68 :
69 : void
70 63 : MFEMCutTransitionSubMesh::labelMesh(mfem::ParMesh & parent_mesh)
71 : {
72 63 : int mpi_comm_rank = getMFEMProblem().getProblemData().myid;
73 63 : int mpi_comm_size = getMFEMProblem().getProblemData().num_procs;
74 :
75 : // First determine face normal based on the first boundary element found on the cut
76 : // to use when determining orientation relative to the cut
77 63 : const mfem::Array<int> & parent_cut_element_id_map = _cut_submesh->GetParentElementIDMap();
78 63 : int rank_with_submesh = -1;
79 63 : if (parent_cut_element_id_map.Size() > 0)
80 : {
81 55 : int reference_face = parent_cut_element_id_map[0];
82 55 : _cut_normal = findFaceNormal(parent_mesh, parent_mesh.GetBdrElementFaceIndex(reference_face));
83 55 : rank_with_submesh = mpi_comm_rank;
84 : }
85 63 : MPI_Allreduce(MPI_IN_PLACE, &rank_with_submesh, 1, MPI_INT, MPI_MAX, getMFEMProblem().getComm());
86 63 : MPI_Bcast(_cut_normal.GetData(),
87 : _cut_normal.Size(),
88 : MFEM_MPI_REAL_T,
89 : rank_with_submesh,
90 : getMFEMProblem().getComm());
91 : // Iterate over all vertices on cut, find elements with those vertices,
92 : // and declare them transition elements if they are on the +ve side of the cut
93 63 : mfem::Array<int> transition_els;
94 63 : std::vector<HYPRE_BigInt> global_cut_vert_ids;
95 63 : mfem::Array<HYPRE_BigInt> gi;
96 63 : parent_mesh.GetGlobalVertexIndices(gi);
97 63 : std::unique_ptr<mfem::Table> vert_to_elem(parent_mesh.GetVertexToElementTable());
98 63 : const mfem::Array<int> & cut_to_parent_vertex_id_map = _cut_submesh->GetParentVertexIDMap();
99 4857 : for (int i = 0; i < _cut_submesh->GetNV(); ++i)
100 : {
101 4794 : int cut_vert = cut_to_parent_vertex_id_map[i];
102 4794 : global_cut_vert_ids.push_back(gi[cut_vert]);
103 4794 : int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
104 4794 : const int * els_adj_to_cut = vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
105 116858 : for (int i = 0; i < ne; i++)
106 : {
107 112064 : const int el_adj_to_cut = els_adj_to_cut[i];
108 204210 : if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
109 92146 : isPositiveSideOfCut(el_adj_to_cut, cut_vert, parent_mesh))
110 46134 : transition_els.Append(el_adj_to_cut);
111 : }
112 : }
113 :
114 : // share cut verts coords or global ids across all procs
115 63 : int n_cut_vertices = global_cut_vert_ids.size();
116 63 : std::vector<int> cut_vert_sizes(mpi_comm_size);
117 63 : MPI_Allgather(
118 : &n_cut_vertices, 1, MPI_INT, cut_vert_sizes.data(), 1, MPI_INT, getMFEMProblem().getComm());
119 : // Make an offset array and total the sizes.
120 63 : std::vector<int> n_vert_offset(mpi_comm_size);
121 63 : std::exclusive_scan(cut_vert_sizes.begin(), cut_vert_sizes.end(), n_vert_offset.begin(), 0);
122 63 : int global_n_cut_vertices = std::accumulate(cut_vert_sizes.begin(), cut_vert_sizes.end(), 0);
123 :
124 : // Gather the queries to all ranks.
125 63 : std::vector<HYPRE_BigInt> all_cut_verts(global_n_cut_vertices);
126 63 : MPI_Allgatherv(global_cut_vert_ids.data(),
127 : n_cut_vertices,
128 : HYPRE_MPI_BIG_INT,
129 : all_cut_verts.data(),
130 : cut_vert_sizes.data(),
131 : n_vert_offset.data(),
132 : HYPRE_MPI_BIG_INT,
133 : getMFEMProblem().getComm());
134 :
135 : // Detect shared vertices and add corresponding elements
136 103 : for (const auto g : make_range(1, parent_mesh.GetNGroups()))
137 22896 : for (const auto gv : make_range(parent_mesh.GroupNVertices(g)))
138 : {
139 : // all elements touching this shared vertex should be updated
140 22856 : int cut_vert = parent_mesh.GroupVertex(g, gv);
141 2447392 : for (std::size_t i = 0; i < all_cut_verts.size(); i += 1)
142 2424536 : if (gi[cut_vert] == all_cut_verts[i]) // check if shared vertex is on the cut plane
143 : {
144 608 : int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
145 : const int * els_adj_to_cut =
146 608 : vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
147 8528 : for (int i = 0; i < ne; i++)
148 : {
149 7920 : const int el_adj_to_cut = els_adj_to_cut[i];
150 15068 : if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
151 7148 : isPositiveSideOfCut(el_adj_to_cut, cut_vert, parent_mesh))
152 3528 : transition_els.Append(el_adj_to_cut);
153 : }
154 : }
155 : }
156 :
157 63 : transition_els.Sort();
158 63 : transition_els.Unique();
159 :
160 63 : setAttributes(parent_mesh, transition_els);
161 63 : }
162 :
163 : mfem::Vector
164 55 : MFEMCutTransitionSubMesh::findFaceNormal(const mfem::ParMesh & mesh, const int & face)
165 : {
166 55 : if (mesh.SpaceDimension() != 3)
167 0 : mooseError("MFEMCutTransitionSubMesh only works in 3-dimensional meshes!");
168 55 : mfem::Vector normal;
169 55 : mfem::Array<int> face_verts;
170 55 : std::vector<mfem::Vector> v;
171 55 : mesh.GetFaceVertices(face, face_verts);
172 :
173 : // First we get the coordinates of 3 vertices on the face
174 220 : for (auto vtx : face_verts)
175 : {
176 165 : mfem::Vector vtx_coords(3);
177 660 : for (int j = 0; j < 3; ++j)
178 495 : vtx_coords[j] = mesh.GetVertex(vtx)[j];
179 165 : v.push_back(vtx_coords);
180 165 : }
181 :
182 : // Now we find the unit vector normal to the face
183 55 : v[0] -= v[1];
184 55 : v[1] -= v[2];
185 55 : v[0].cross3D(v[1], normal);
186 55 : normal /= normal.Norml2();
187 110 : return normal;
188 55 : }
189 :
190 : bool
191 99294 : MFEMCutTransitionSubMesh::isPositiveSideOfCut(const int & el,
192 : const int & el_vertex_on_cut,
193 : mfem::ParMesh & parent_mesh)
194 : {
195 99294 : const int sdim = parent_mesh.SpaceDimension();
196 99294 : mfem::Vector el_center(3);
197 99294 : parent_mesh.GetElementCenter(el, el_center);
198 99294 : mfem::Vector vertex_coords(parent_mesh.GetVertex(el_vertex_on_cut), sdim);
199 99294 : mfem::Vector relative_center(sdim);
200 397176 : for (int j = 0; j < sdim; j++)
201 297882 : relative_center[j] = el_center[j] - vertex_coords[j];
202 198588 : return _cut_normal * relative_center > 0;
203 99294 : }
204 :
205 : void
206 63 : MFEMCutTransitionSubMesh::setAttributes(mfem::ParMesh & parent_mesh,
207 : mfem::Array<int> & transition_els)
208 : {
209 : // Generate a set of local new attributes for transition region elements
210 63 : const int old_max_attr = parent_mesh.attributes.Max();
211 63 : mfem::Array<int> new_attrs(old_max_attr);
212 63 : new_attrs = -1;
213 23354 : for (auto const & transition_el : transition_els)
214 : {
215 23291 : const int old_attr = parent_mesh.GetAttribute(transition_el);
216 23291 : new_attrs[old_attr - 1] = old_max_attr + old_attr;
217 : // Set the new attribute IDs for transition region elements
218 23291 : parent_mesh.SetAttribute(transition_el, new_attrs[old_attr - 1]);
219 : }
220 :
221 : // Distribute local attribute IDs to other MPI ranks to create set of new
222 : // global attribute IDs for the transition region.
223 63 : MPI_Allreduce(
224 : MPI_IN_PLACE, new_attrs, old_max_attr, MPI_INT, MPI_MAX, getMFEMProblem().getComm());
225 :
226 63 : mfem::AttributeSets & attr_sets = parent_mesh.attribute_sets;
227 : // Create attribute set labelling the newly created transition region on one side of the cut
228 63 : attr_sets.CreateAttributeSet(_transition_subdomain);
229 : // Create attribute set labelling the entire closed geometry
230 63 : attr_sets.SetAttributeSet(_closed_subdomain, getSubdomainAttributes());
231 : // Add the new domain attributes to new attribute sets
232 63 : const std::set<std::string> attr_set_names = attr_sets.GetAttributeSetNames();
233 252 : for (int old_attr = 1; old_attr <= old_max_attr; ++old_attr)
234 : {
235 189 : int new_attr = new_attrs[old_attr - 1];
236 : // Add attributes only if they're transition region attributes
237 189 : if (new_attr != -1)
238 : {
239 101 : attr_sets.AddToAttributeSet(_transition_subdomain, new_attr);
240 606 : for (auto const & attr_set_name : attr_set_names)
241 : {
242 505 : const mfem::Array<int> & attr_set = attr_sets.GetAttributeSet(attr_set_name);
243 : // If attribute set has the old attribute of the transition region, add the new one
244 505 : if (attr_set.Find(old_attr) != -1)
245 202 : attr_sets.AddToAttributeSet(attr_set_name, new_attr);
246 : }
247 : }
248 : }
249 :
250 63 : parent_mesh.SetAttributes();
251 63 : }
252 :
253 : bool
254 119984 : MFEMCutTransitionSubMesh::isInDomain(const int & element,
255 : const mfem::Array<int> & subdomains,
256 : const mfem::ParMesh & mesh)
257 : {
258 : // element<0 for ghost elements
259 119984 : if (element < 0)
260 0 : return true;
261 :
262 166776 : for (const auto & subdomain : subdomains)
263 146086 : if (mesh.GetAttribute(element) == subdomain)
264 99294 : return true;
265 20690 : return false;
266 : }
267 :
268 : #endif
|
