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 "surfacecurrentoperator.hpp"
5 :
6 : #include "fem/coefficient.hpp"
7 : #include "utils/communication.hpp"
8 : #include "utils/geodata.hpp"
9 : #include "utils/iodata.hpp"
10 :
11 : namespace palace
12 : {
13 :
14 0 : SurfaceCurrentData::SurfaceCurrentData(const config::SurfaceCurrentData &data,
15 : const mfem::ParMesh &mesh)
16 : {
17 : // Construct the source elements allowing for a possible multielement surface current
18 : // sources.
19 0 : for (const auto &elem : data.elements)
20 : {
21 : mfem::Array<int> attr_list;
22 0 : attr_list.Append(elem.attributes.data(), elem.attributes.size());
23 0 : switch (elem.coordinate_system)
24 : {
25 0 : case CoordinateSystem::CYLINDRICAL:
26 0 : elems.push_back(
27 0 : std::make_unique<CoaxialElementData>(elem.direction, attr_list, mesh));
28 0 : break;
29 0 : case CoordinateSystem::CARTESIAN:
30 0 : elems.push_back(
31 0 : std::make_unique<UniformElementData>(elem.direction, attr_list, mesh));
32 0 : break;
33 : }
34 : }
35 0 : }
36 :
37 0 : double SurfaceCurrentData::GetExcitationCurrent() const
38 : {
39 : // Ideal unit current source for each index.
40 0 : return 1.0;
41 : }
42 :
43 20 : SurfaceCurrentOperator::SurfaceCurrentOperator(const IoData &iodata,
44 : const mfem::ParMesh &mesh)
45 : {
46 : // Set up surface current source boundaries.
47 20 : SetUpBoundaryProperties(iodata, mesh);
48 20 : PrintBoundaryInfo(iodata, mesh);
49 20 : }
50 :
51 20 : void SurfaceCurrentOperator::SetUpBoundaryProperties(const IoData &iodata,
52 : const mfem::ParMesh &mesh)
53 : {
54 : // Check that surface current boundary attributes have been specified correctly.
55 20 : if (!iodata.boundaries.current.empty())
56 : {
57 0 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
58 : mfem::Array<int> bdr_attr_marker(bdr_attr_max), source_marker(bdr_attr_max);
59 : bdr_attr_marker = 0;
60 : source_marker = 0;
61 0 : for (auto attr : mesh.bdr_attributes)
62 : {
63 0 : bdr_attr_marker[attr - 1] = 1;
64 : }
65 0 : for (const auto &[idx, data] : iodata.boundaries.current)
66 : {
67 0 : for (const auto &elem : data.elements)
68 : {
69 0 : for (auto attr : elem.attributes)
70 : {
71 0 : MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
72 : "Surface current boundary attribute tags must be non-negative and "
73 : "correspond to boundaries in the mesh!");
74 0 : MFEM_VERIFY(bdr_attr_marker[attr - 1],
75 : "Unknown surface current boundary attribute " << attr << "!");
76 0 : MFEM_VERIFY(
77 : !source_marker[attr - 1],
78 : "Boundary attribute is assigned to more than one surface current source!");
79 0 : source_marker[attr - 1] = 1;
80 : }
81 : }
82 : }
83 : }
84 :
85 : // Set up surface current data structures.
86 20 : for (const auto &[idx, data] : iodata.boundaries.current)
87 : {
88 0 : sources.try_emplace(idx, data, mesh);
89 : }
90 20 : }
91 :
92 20 : void SurfaceCurrentOperator::PrintBoundaryInfo(const IoData &iodata,
93 : const mfem::ParMesh &mesh)
94 : {
95 20 : if (sources.empty())
96 : {
97 : return;
98 : }
99 0 : Mpi::Print("\nConfiguring surface current excitation source term at attributes:\n");
100 0 : for (const auto &[idx, data] : sources)
101 : {
102 0 : for (const auto &elem : data.elems)
103 : {
104 0 : for (auto attr : elem->GetAttrList())
105 : {
106 0 : Mpi::Print(" {:d}: Index = {:d}, n = ({:+.1f})\n", attr, idx,
107 0 : fmt::join(mesh::GetSurfaceNormal(mesh, attr), ","));
108 : }
109 : }
110 : }
111 : }
112 :
113 0 : const SurfaceCurrentData &SurfaceCurrentOperator::GetSource(int idx) const
114 : {
115 : auto it = sources.find(idx);
116 0 : MFEM_VERIFY(it != sources.end(), "Unknown current source index requested!");
117 0 : return it->second;
118 : }
119 :
120 20 : mfem::Array<int> SurfaceCurrentOperator::GetAttrList() const
121 : {
122 : mfem::Array<int> attr_list;
123 20 : for (const auto &[idx, data] : sources)
124 : {
125 0 : for (const auto &elem : data.elems)
126 : {
127 : attr_list.Append(elem->GetAttrList());
128 : }
129 : }
130 20 : return attr_list;
131 : }
132 :
133 0 : void SurfaceCurrentOperator::AddExcitationBdrCoefficients(SumVectorCoefficient &fb)
134 : {
135 : // Construct the RHS source term for surface current boundaries, which looks like
136 : // -iω J_inc for a surface current boundary. The chosen surface current J_inc corresponds
137 : // to a unit current excitation. Note: The real RHS returned here does not yet have the
138 : // factor of (iω) included, so works for time domain simulations requiring RHS -J_inc
139 : // (t).
140 0 : for (const auto &[idx, data] : sources)
141 : {
142 0 : AddExcitationBdrCoefficients(data, fb);
143 : }
144 0 : }
145 :
146 0 : void SurfaceCurrentOperator::AddExcitationBdrCoefficients(int idx, SumVectorCoefficient &fb)
147 : {
148 : // Construct the RHS source term for a single surface current boundary index.
149 0 : AddExcitationBdrCoefficients(GetSource(idx), fb);
150 0 : }
151 :
152 0 : void SurfaceCurrentOperator::AddExcitationBdrCoefficients(const SurfaceCurrentData &data,
153 : SumVectorCoefficient &fb)
154 : {
155 : // Add excited boundaries to the linear form, with a unit current distributed across
156 : // all elements of the current source in parallel.
157 0 : for (const auto &elem : data.elems)
158 : {
159 0 : const double Jinc = 1.0 / (elem->GetGeometryWidth() * data.elems.size());
160 0 : fb.AddCoefficient(elem->GetModeCoefficient(-Jinc));
161 : }
162 0 : }
163 :
164 : } // namespace palace
|
