Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_UTILS_GEODATA_HPP
5 : #define PALACE_UTILS_GEODATA_HPP
6 :
7 : #include <array>
8 : #include <cmath>
9 : #include <memory>
10 : #include <vector>
11 : #include <mfem.hpp>
12 :
13 : namespace palace
14 : {
15 :
16 : class IoData;
17 :
18 : namespace mesh
19 : {
20 :
21 : //
22 : // Functions for mesh related functionality.
23 : //
24 :
25 : // Read and partition a serial mesh from file, returning a pointer to the new parallel mesh
26 : // object, which should be destroyed by the user.
27 : std::unique_ptr<mfem::ParMesh> ReadMesh(const IoData &iodata, MPI_Comm comm);
28 :
29 : // Refine the provided mesh according to the data in the input file. If levels of refinement
30 : // are requested, the refined meshes are stored in order of increased refinement. Ownership
31 : // of the initial coarse mesh is inherited by the fine meshes and it should not be deleted.
32 : // The fine mesh hierarchy is owned by the user.
33 : void RefineMesh(const IoData &iodata, std::vector<std::unique_ptr<mfem::ParMesh>> &mesh);
34 :
35 : // Dimensionalize a mesh for use in exporting a mesh. Scales vertices and nodes by L.
36 : void DimensionalizeMesh(mfem::Mesh &mesh, double L);
37 :
38 : // Nondimensionalize a mesh for use in the solver. Scales vertices and nodes by 1/L.
39 : void NondimensionalizeMesh(mfem::Mesh &mesh, double L);
40 :
41 : // Struct containing flags for the (global) mesh element types.
42 : struct ElementTypeInfo
43 : {
44 : bool has_simplices;
45 : bool has_hexahedra;
46 : bool has_prisms;
47 : bool has_pyramids;
48 : std::vector<mfem::Geometry::Type> GetGeomTypes() const;
49 : };
50 :
51 : // Simplified helper for describing the element types in a (Par)Mesh.
52 : ElementTypeInfo CheckElements(const mfem::Mesh &mesh);
53 :
54 : // Check if a tetrahedral (Par)Mesh is ready for local refinement.
55 : bool CheckRefinementFlags(const mfem::Mesh &mesh);
56 :
57 : // Helper function to convert a set of attribute numbers to a marker array. The marker array
58 : // will be of size max_attr and it will contain only zeroes and ones. Ones indicate which
59 : // attribute numbers are present in the list array. In the special case when list has a
60 : // single entry equal to -1 the marker array will contain all ones.
61 : void AttrToMarker(int max_attr, const int *attr_list, int attr_list_size,
62 : mfem::Array<int> &marker, bool skip_invalid = false);
63 :
64 : template <typename T>
65 : inline void AttrToMarker(int max_attr, const T &attr_list, mfem::Array<int> &marker,
66 : bool skip_invalid = false)
67 : {
68 : const auto size = std::distance(attr_list.begin(), attr_list.end());
69 473 : AttrToMarker(max_attr, (size > 0) ? &attr_list[0] : nullptr, size, marker, skip_invalid);
70 317 : }
71 :
72 : template <typename T>
73 311 : inline mfem::Array<int> AttrToMarker(int max_attr, const T &attr_list,
74 : bool skip_invalid = false)
75 : {
76 : mfem::Array<int> marker;
77 : AttrToMarker(max_attr, attr_list, marker, skip_invalid);
78 311 : return marker;
79 : }
80 :
81 : // Helper function to construct the axis-aligned bounding box for all elements with the
82 : // given attribute.
83 : void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
84 : bool bdr, mfem::Vector &min, mfem::Vector &max);
85 :
86 : inline void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr,
87 : mfem::Vector &min, mfem::Vector &max)
88 : {
89 : mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
90 : marker = 0;
91 : marker[attr - 1] = 1;
92 : GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max);
93 : }
94 :
95 73 : inline void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, mfem::Vector &min,
96 : mfem::Vector &max)
97 : {
98 73 : mfem::Array<int> marker(mesh.attributes.Max());
99 : marker = 1;
100 73 : GetAxisAlignedBoundingBox(mesh, marker, false, min, max);
101 73 : }
102 :
103 : // Struct describing a bounding box in terms of the center and face normals. The normals
104 : // specify the direction from the center of the box.
105 : struct BoundingBox
106 : {
107 : // The central point of the bounding box.
108 : std::array<double, 3> center;
109 :
110 : // Vectors from center to the midpoint of each face.
111 : std::array<std::array<double, 3>, 3> axes;
112 :
113 : // Whether or not this bounding box is two dimensional.
114 : bool planar;
115 :
116 : // Compute the area of the bounding box spanned by the first two normals.
117 : double Area() const;
118 :
119 : // Compute the volume of the 3D bounding box. Returns zero if planar.
120 : double Volume() const;
121 :
122 : // Compute the normalized axes of the bounding box.
123 : std::array<std::array<double, 3>, 3> Normals() const;
124 :
125 : // Compute the lengths along each axis.
126 : std::array<double, 3> Lengths() const;
127 :
128 : // Compute the deviations in degrees of a vector from each of the axis directions. Angles
129 : // are returned in the interval [0, 180].
130 : std::array<double, 3> Deviations(const std::array<double, 3> &direction) const;
131 : };
132 :
133 : // Helper functions for computing bounding boxes from a mesh and markers. These do not need
134 : // to be axis-aligned. Note: This function only returns a minimum oriented bounding box for
135 : // points whose convex hull exactly forms a rectangle or rectangular prism, implementing a
136 : // vastly simplified version of QuickHull for this case. For other shapes, the result is
137 : // less predictable, and may not even form a bounding box of the sampled point cloud.
138 : BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
139 : bool bdr);
140 :
141 : inline BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr)
142 : {
143 : mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
144 : marker = 0;
145 : marker[attr - 1] = 1;
146 : return GetBoundingBox(mesh, marker, bdr);
147 : }
148 :
149 : // Given a mesh and a marker, compute the bounding circle/sphere of the marked elements. In
150 : // this case the normals of the bounding box object are arbitrary, and the Area and Volume
151 : // members should not be used, but the Lengths function returns the ball diameter. This
152 : // function implements Welzl's algorithm.
153 : BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
154 : bool bdr);
155 :
156 : inline BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, int attr, bool bdr)
157 : {
158 : mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
159 : marker = 0;
160 : marker[attr - 1] = 1;
161 : return GetBoundingBall(mesh, marker, bdr);
162 : }
163 :
164 : // Helper function for computing the direction aligned length of a marked group.
165 : double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
166 : bool bdr, const std::array<double, 3> &dir);
167 :
168 : inline double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr,
169 : const std::array<double, 3> &dir)
170 : {
171 : mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
172 : marker = 0;
173 : marker[attr - 1] = 1;
174 : return GetProjectedLength(mesh, marker, bdr, dir);
175 : }
176 :
177 : // Helper function for computing the closest distance of a marked group to a given point,
178 : // by brute force searching over the entire point set. Optionally compute the furthest
179 : // distance instead of the closest.
180 : double GetDistanceFromPoint(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
181 : bool bdr, const std::array<double, 3> &origin,
182 : bool max = false);
183 :
184 : inline double GetDistanceFromPoint(const mfem::ParMesh &mesh, int attr, bool bdr,
185 : const std::array<double, 3> &dir, bool max = false)
186 : {
187 : mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
188 : marker = 0;
189 : marker[attr - 1] = 1;
190 : return GetDistanceFromPoint(mesh, marker, bdr, dir, max);
191 : }
192 :
193 : // Helper function to compute the average surface normal for all elements with the given
194 : // attributes.
195 : mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
196 : bool average = true);
197 :
198 78 : inline mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, int attr,
199 : bool average = true)
200 : {
201 : const bool bdr = (mesh.Dimension() == mesh.SpaceDimension());
202 78 : mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
203 : marker = 0;
204 78 : marker[attr - 1] = 1;
205 156 : return GetSurfaceNormal(mesh, marker, average);
206 : }
207 :
208 0 : inline mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, bool average = true)
209 : {
210 : const bool bdr = (mesh.Dimension() == mesh.SpaceDimension());
211 0 : const auto &attributes = bdr ? mesh.bdr_attributes : mesh.attributes;
212 0 : return GetSurfaceNormal(mesh, AttrToMarker(attributes.Max(), attributes), average);
213 : }
214 :
215 : // Helper functions to compute the volume or area for all domain or boundary elements with
216 : // the given attributes.
217 : double GetSurfaceArea(const mfem::ParMesh &mesh, const mfem::Array<int> &marker);
218 :
219 : inline double GetSurfaceArea(const mfem::ParMesh &mesh, int attr)
220 : {
221 : mfem::Array<int> marker(mesh.bdr_attributes.Max());
222 : marker = 0;
223 : marker[attr - 1] = 1;
224 : return GetSurfaceArea(mesh, marker);
225 : }
226 :
227 : double GetVolume(const mfem::ParMesh &mesh, const mfem::Array<int> &marker);
228 :
229 : inline double GetVolume(const mfem::ParMesh &mesh, int attr)
230 : {
231 : mfem::Array<int> marker(mesh.attributes.Max());
232 : marker = 0;
233 : marker[attr - 1] = 1;
234 : return GetVolume(mesh, marker);
235 : }
236 :
237 : // Helper function responsible for rebalancing the mesh, and optionally writing meshes from
238 : // the intermediate stages to disk. Returns the imbalance ratio before rebalancing.
239 : double RebalanceMesh(const IoData &iodata, std::unique_ptr<mfem::ParMesh> &mesh);
240 :
241 : // Helper for creating a hexahedral mesh from a tetrahedral mesh.
242 : mfem::Mesh MeshTetToHex(const mfem::Mesh &orig_mesh);
243 :
244 : } // namespace mesh
245 :
246 : } // namespace palace
247 :
248 : #endif // PALACE_UTILS_GEODATA_HPP
|
