Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "lumpedelement.hpp"
5 :
6 : #include "fem/coefficient.hpp"
7 : #include "fem/integrator.hpp"
8 : #include "linalg/vector.hpp"
9 : #include "utils/communication.hpp"
10 : #include "utils/geodata.hpp"
11 :
12 : namespace palace
13 : {
14 :
15 30 : UniformElementData::UniformElementData(const std::array<double, 3> &input_dir,
16 : const mfem::Array<int> &attr_list,
17 30 : const mfem::ParMesh &mesh)
18 30 : : LumpedElementData(attr_list)
19 : {
20 30 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
21 30 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
22 30 : auto bounding_box = mesh::GetBoundingBox(mesh, attr_marker, true);
23 :
24 : // Check the user specified direction aligns with an axis direction.
25 30 : constexpr double angle_warning_deg = 0.1;
26 : constexpr double angle_error_deg = 1.0;
27 30 : auto lengths = bounding_box.Lengths();
28 30 : auto deviations_deg = bounding_box.Deviations(input_dir);
29 30 : if (std::none_of(deviations_deg.begin(), deviations_deg.end(),
30 : [](double x) { return x < angle_warning_deg; }))
31 : {
32 0 : auto normals = bounding_box.Normals();
33 0 : Mpi::Warning("User specified direction {} does not align with either bounding box "
34 : "axis up to {:.3e} degrees!\n"
35 : "Axis 1: {} ({:.3e} degrees)\nAxis 2: {} ({:.3e} degrees)\nAxis 3: "
36 : "{} ({:.3e} degrees)!\n",
37 : input_dir, angle_warning_deg, normals[0], deviations_deg[0], normals[1],
38 : deviations_deg[1], normals[2], deviations_deg[2]);
39 : }
40 30 : if (std::none_of(deviations_deg.begin(), deviations_deg.end(),
41 : [](double x) { return x < angle_error_deg; }))
42 : {
43 : Mpi::Barrier(mesh.GetComm());
44 0 : MFEM_ABORT("Specified direction does not align sufficiently with bounding box axes ("
45 : << deviations_deg[0] << ", " << deviations_deg[1] << ", "
46 : << deviations_deg[2] << " vs. tolerance " << angle_error_deg << ")!");
47 : }
48 30 : direction.SetSize(input_dir.size());
49 : std::copy(input_dir.begin(), input_dir.end(), direction.begin());
50 30 : direction /= direction.Norml2();
51 :
52 : // Compute the length from the most aligned normal direction.
53 : constexpr double rel_tol = 1.0e-6;
54 : auto l_component =
55 : std::distance(deviations_deg.begin(),
56 : std::min_element(deviations_deg.begin(), deviations_deg.end()));
57 30 : l = lengths[l_component];
58 30 : MFEM_VERIFY(std::abs(l - mesh::GetProjectedLength(mesh, attr_marker, true, input_dir)) <
59 : rel_tol * l,
60 : "Bounding box discovered length ("
61 : << l << ") should match projected length ("
62 : << mesh::GetProjectedLength(mesh, attr_marker, true, input_dir) << "!");
63 :
64 : // Compute the width as area / length. This allows the lumped element to be non-planar,
65 : // and generalizes nicely to the case for an infinitely thin rectangular lumped element
66 : // with elements on both sides (for which the width computed from the bounding box would
67 : // be incorrect by a factor of 2).
68 30 : double area = mesh::GetSurfaceArea(mesh, attr_marker);
69 30 : MFEM_VERIFY(area > 0.0, "Uniform lumped element has zero area!");
70 30 : w = area / l;
71 30 : }
72 :
73 : std::unique_ptr<mfem::VectorCoefficient>
74 12 : UniformElementData::GetModeCoefficient(double coeff) const
75 : {
76 12 : mfem::Vector source = direction;
77 12 : source *= coeff;
78 12 : return std::make_unique<RestrictedVectorCoefficient<mfem::VectorConstantCoefficient>>(
79 24 : attr_list, source);
80 : }
81 :
82 0 : CoaxialElementData::CoaxialElementData(const std::array<double, 3> &input_dir,
83 : const mfem::Array<int> &attr_list,
84 0 : const mfem::ParMesh &mesh)
85 0 : : LumpedElementData(attr_list)
86 : {
87 0 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
88 0 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
89 0 : auto bounding_ball = mesh::GetBoundingBall(mesh, attr_marker, true);
90 0 : MFEM_VERIFY(bounding_ball.planar,
91 : "Boundary elements must be coplanar to define a coaxial lumped element!");
92 :
93 : // Direction of the excitation as +/-r̂.
94 0 : direction = (input_dir[0] > 0 ? +1 : -1);
95 0 : origin.SetSize(bounding_ball.center.size());
96 : std::copy(bounding_ball.center.begin(), bounding_ball.center.end(), origin.begin());
97 :
98 : // Get outer and inner radius of the annulus, assuming full 2π circumference.
99 0 : r_outer = 0.5 * bounding_ball.Lengths()[0];
100 0 : r_inner = mesh::GetDistanceFromPoint(mesh, attr_marker, true, bounding_ball.center);
101 0 : MFEM_VERIFY(r_inner > 0.0,
102 : "Coaxial element annulus has should have positive inner radius!");
103 0 : MFEM_VERIFY(r_outer > r_inner, "Coaxial element annulus has unexpected outer radius "
104 : << r_outer << " <= inner radius " << r_inner << "!");
105 0 : }
106 :
107 : std::unique_ptr<mfem::VectorCoefficient>
108 0 : CoaxialElementData::GetModeCoefficient(double coeff) const
109 : {
110 0 : coeff *= direction;
111 0 : mfem::Vector x0(origin);
112 0 : auto Source = [coeff, x0](const mfem::Vector &x, mfem::Vector &f) -> void
113 : {
114 0 : f = x;
115 0 : f -= x0;
116 0 : double oor = 1.0 / f.Norml2();
117 0 : f *= coeff * oor * oor;
118 0 : };
119 0 : return std::make_unique<RestrictedVectorCoefficient<mfem::VectorFunctionCoefficient>>(
120 0 : attr_list, x0.Size(), Source);
121 : }
122 :
123 : } // namespace palace
|
